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ABSTRACT 



The Naval Weapons Support Center is planning to implement a bonus system to 
improve the reliability of pyrotechnic devices. The measure of effectiveness that they 
wish to use to determine how to award bonuses is the reliability of pyrotechnic devices. 
The data available to estimate this reliability is based on the current sampling inspection 
plan in which devices are tested in different environments. The models which include 
both dependence and independence assumptions between the outcomes of these tests are 
implemented and estimates of overall reliability along with 95 % lower confidence bound 
are obtained. The 95 % lower confidence bounds are found by bootstrapping. Using 
these estimates, models for making the decision to award bonuses are discussed and 
studied using Monte Carlo simulation . 
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I. INTRODUCTION 



A. BACKGROUND 

Nations spend a lot of money to establish a strong defense network. It is essential 
that nations buy reliable weapons and ammunition from the contractors. To ensure re- 
liability, contracts must include a lot acceptance sampling plan that specifies the minimal 
acceptable quality. 

The Naval Weapons Support Center purchases pyrotechnic devices. Unless other- 
wise specified in the contract, the supplier is responsible to see that his devices meet all 
inspection requirements as specified. The inspection requirements are particular to 
characteristics of each type of device and are specified on the reference drawings and 
supplemental quality assurance provisions of the contract. Any testing that needs to be 
done on these devices is explained in this contract. 

When nations buy weapons and ammunition from the same contractor, they would 
like the quality to improve over time. As contracts are now written, contractors need 
only to satisfy the requirements of the sampling inspection plan for lot acceptance. 
Under such contracts, contractors have no incentive to improve the quality of items they 
provide. For this reason, to improve quality. The Naval Weapons Support Center has 
decided to implement a bonus system. The contractor will be awarded a bonus if the 
result of the sampling inspection exceeds the minimum requirements for lot acceptance. 

The Naval Weapons Support Center will begin to implement a bonus system for 
pyrotechnic devices in FY91. The data available to make the decision whether to award 
a bonus is based on the current sampling inspection plan. This plan is a series of de- 
structive tests in different environments. The purpose of this thesis is to provide the 
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Naval Weapons Support Center with guidance for implementing a bonus system for 
pyrotechnic devices. Once implemented, the pyrotechnic bonus system will serve as a 
prototype for bonus systems for other devices. 

B. PYROTECHNIC DEVICE RELIABILITY 

A pyrotechnic device is a chemical and grenade ammunition. There are three cate- 
gories of pyrotechnic devices. The three categories are: Aerial display, Surface display 
and Grenades [Ref 1: pp. I-2j. Samples of the Pyrotechnic devices are exposed to 
various environments and then activated. The criteria for successful activation depends 
on the type of device : 

1. Aerial display (Ground signals, flares, airburst simulators, sub signals, signal kits, 

etc.) 

"Successful activation means that the item will, after simulating user environment, 
successfully " 

• launch, 

• have proper separation / signal ignition, 

• reach desired altitude at correct angle, 

• have proper parachute deployment, 

• have proper display color, 

• have proper display time, 

• have no subsequent interference with next item. 

2. Surface (Ground or Water) display (flare, hand-held signals, smoke / illume 
grenades, simulators, etc.) 

"Successful activation means that the item will, after simulating user environment, 
successfully ..." 
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• have signal ignition with proper display, 

• have proper display, 

• have proper display time, 

3. Grenades (fragmentation, defensive, white phosphors, etc.) 

"Successful activation means that the item will, after simulating user environment, 
successfully 

• function (high order detonation following delay), 

• have proper dissipation of payload, 

• have completely consumed payload. 

C. PYROTECHNIC DEVICE RELIABILITY PROBLEM 

Samples from any large lot of pyrotechnic devices submitted by a manufacturer must 
activate after exposure to different environments. These are; 

1. Manufacturer Environment, 

2. Temperature and Humidity Environment, 

3. Vibration Environment, 

4. Altitude Environment. 

All items tested are subjected to the manufacturer environment. However items are only 
subjected to one of the three remaining environments : Temperature and Humidity, 
Vibration or Altitude. 

The sampling plan consists of using four distinct samples from a lot that can be as- 
sumed (approximately) statistically independent. The items tested in each sample are 
also assumed to be independent. According to the sampling plan; 

• 20 items are subjected to the Manufacturer Test, 
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• 20 items are subjected to both the Temperature and Humidity Test and Man- 
ufacturer Test, 

• 32 items are subjected to both the Vibration Test and Manufacturer Test, 

• 20 items are subjected to both Altitude Test and Manufacturer Test. 

A total of 92 items are tested. 

Acceptance criteria for each test are : 

1. Manufacturer Test : Of the 20 items; if no more than 1 fails to activate, the lot 
passes. 

2. Joint Temperature and Humidity and Manufacturer Test : Of the 20 items; if no 
more than 1 fails to activate, the lot passes. 

3. Joint Vibration and Manufacturer Test : Of the 32 items; if no more than 2 fails 
to activate, the lot passes. 

4. Joint Altitude and Manufacturer Test : Of the 20 items; if no more than 1 fails to 
activate, the lot passes. 

The number of failures for these tests will be summarized by the vector : 

( FOM, FOTH, FOV, FOA ) (1.2) 



where 

• FOM represents the number of failures after the manufacturer test, 

• FOTH represents the number of failures after both the temperature and humidity 
and manufacturer tests, 

• FOV represents the number of failures after both the vibration and manufacturer 
tests, 

• FOA represents the number of failures after the altitude and manufacturer tests. 
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For example. 



( 1 , 1 , 2 , 1 ). ( 1 . 2 ) 

represents, 1 failure in manufacturer test, 1 failure in the joint temperature and humidity 
and manufacturer test, 2 failures in the joint vibration and manufacturer test, and 1 
failure in the joint altitude and manufacturer test. It is also the maximum number of 
failures in each test that still leads to lot acceptance. 

The marginal distributions of the number of failures in each of the above tests are 
modeled by binomial distributions. There are 24 possible realizations of the sampling 
inspection; ranging from the best case with no failures to activate, to the worst case with 
1, 1, 2 and 1 devices failing to activate in these tests manufacturer, temperature and 
humidity, vibration and altitude respectively. These cases are tabulated in Table 1. 

To award bonuses we need one measure of effectiveness for pyrotechnic devices that 
can be estimated from the available data. Ideally, this measure is the reliability of the 
device. However, because each of the tests are destructive, there is no one natural defi- 
nition of reliability for these devices. To be on the conservative side, we define the reli- 
ability of a device to be the probability that the device will activate after exposure to all 
of the environments. It is not easy to estimate this reliability from the sampling plan 
data. This data is incomplete in the sense that we have limited information about the 
joint probability of activation after exposure to more than one environment. To try to 
compensate for this lack in the data, we will use models for the joint distribution of 
{ FOM, FOTH, FOV, FOA ) that specify particular types of dependence between the 
events that a device activates after exposure to different environments. 
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Table 1. POSSIBLE CASES 



CASE #'S 


FOM 


FOTH 


FOV 


FOA 


1 


0 


0 


0 


0 




0 


0 


0 


1 


3 


0 


0 


1 


0 


4 


0 


1 


0 


0 


5 


1 


0 


0 


0 


6 


0 


0 


1 


1 


7 


0 


1 


0 


1 


8 


0 


1 


1 


0 


9 


1 


0 


0 


1 


10 


1 


0 


1 


0 


11 


1 


1 


0 


0 


12 


0 


0 


2 


0 


13 


0 


1 


1 


1 


14 


1 


0 


1 


1 


15 


1 


1 


0 


1 


16 


1 


1 


1 


0 


17 


0 


0 


2 


1 


18 


0 


1 


2 


0 


19 


1 


0 


2 


0 


20 


1 


1 


1 


1 


21 


0 


1 


2 


1 


22 


1 


0 


2 


1 


23 


1 


1 


2 


0 


24 


1 


1 


2 


1 
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Using these models and based on sampling plan data, estimates of the overall reli- 
ability along with lower confidence bounds are obtained. These will be used to imple- 
ment the bonus system for pyrotechnic devices. We compute the maximum likelihood 
estimator ( MLE ) of the reliability by maximizing the appropriate likelihood; lower 
confidence bounds ( LCB ) are found by bootstrapping. The MLE's are computed under 
both independence and dependence assumptions. The estimation procedures assuming 
independence are described in Chapter II. In Chapter III we incorporate dependence 
by fitting a Log Linear Model to our data. The MLE's from Chapter II and Chapter 
III lead to inappropriate results for this for this problem; thus in Chapter IV we consider 
alternate and ver>' conservative estimates of reliability. Using the estimates of Chapter 
IV, we investigate a sequential scheme for making the decision to award bonuses in 
Chapter V. The results of simulations are presented in Chapter VI. Finally, conclusions 
and recommendations are given in Chapter VII. 
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II. THE ESTIMATION OF THE MAXIMUM LIKELIHOOD ESTIMATOR 



(MLE) OF THE RELIABILITY WITH INDEPENDENCE ASSUMPTION 

A. DEFINITIONS 

We will say that a device survives environment E , if it is still potentially capable of 
activation after exposure to environment E . Let, 

• £, be the device activates after exposure to Manufacturer environment, 

• be the device activates after exposure to Temperature and Humidity environ- 
ment, 

• be the device activates after exposure to Vibration environment, 

• £4 be the device activates after exposure to Altitude environment. 

We define the reliability of device as below, 

£ = p(£, n £2 n ^3 n ^4). (2.1) 

In this formula, R means the probability that a device activates after exposure to four 
environments. We will estimate R for each of the 24 cases which lead to lot acceptance. 

Let Qi = P{Ef) be the probability that device activates after exposure to environ- 
ment i, for i = 1, 2, 3, 4. In the acceptance sampling plan several of the items must 

activate after exposure to a joint manufacturer and another environment. To avoid 
confusion we will denote tests 1 through 4 as the manufacturer test, the joint temper- 
ature humidity and manufacturer test, the joint vibration and manufacturer test and 
joint altitude and manufacturer test respectively. 
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Let 



R, = P{E,) = Q, 



( 2 . 2 ) 



and lei 



/?, = p(£,n £■<•) • 



(2.3) 



Here R, is the probability that a device survives test i for i = 1, 2, 3, 4. The simplest 
model is to assume that £■,, E^, E^, £!, are independent. If we assume that £4 

are independent then 



B. THE LIKELIHOOD EQUATION 

Let 

• be the number of devices that activate after test i, 

• rt, be the number of items given test i. 

Then is binomial with parameters £, and for i = 1, 2, 3, 4. Under the assumption 
of independence the joint likelihood function of observing A’, = jc, , . . . , T 4 = JC 4 is 



Ri = 2 , Qi 



(2.4) 



for i = 2, 3, 4 and the reliability of device is, 



R= Q1Q2 Q3 Qa • 



(2.5) 



4 




( 2 . 6 ) 
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with constraints : 



0 < ^, < 1 , 



0 ^ » 



0 < < R^, 



0 < R4 < Ri. 

Our aim is to maximize this likelihood function subject to the constraints that 
( /?, , , /?j , R,) e S where 

S = {(R^R^^R^R^): 0 ^ R^ < I, 0 < Ri < R^ / = 2, 3, 4 } . 

From the equation (2.6), we see that maximizing L is equivalent to maximizing 



/ = ^{ ( w, In Ri) + (n,- xi) ln( 1 - /?,) }, (2.7) 

/=! 



where the constant multipliers 



for i = 1, 2, 3, 4 have been dropped ( because they 



W 

do not effect the maximization procedure ) and the natural logarithm of L is taken. 

We first show that / is a concave function. To show that / is a concave function, 
we can show - / is a convex function. According to Theorem 3.3.6 [Ref 2: p. 92], 
by looking at its Hessian matrix, we can learn whether function is convex or not. If its 
Hessian matrix is positive semi-definite at each point S then function / is convex. To 
create the Hessian matrix, we must calculate partial derivatives of the function - / , 



dl ^ nj - Xj 

dRi Ri 1 - Ri 



( 2 . 8 ) 
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I 



(2.9) 







( 2 . 10 ) 



Then, the determinant of the Hessian is: 



4 




( 2 . 11 ) 



Clearly we can see that for 0 < ^, < 1 i = 1, 2, 3, 4, | // | is always positive. Because 
- / is continuous, this implies that - / is a convex function on S. As a result of this, 

/ is a concave function. 

We note that with the constraints on the probabilities R, , . . . , Jl, , there does not 
in general exist a closed form solution to MLE. However, with only 24 realizations of 
X, , . . . , X4 of interest, the estimated reliabilities for these 24 cases can be found with 
some rather tedious but straight-forward computations. 

C. COMPUTING THE MAXIMUM LIKELIHOOD ESTIMATORS AND LOWER 
CONFIDENCE BOUNDS 

Because / is concave over the convex set S, if the maximum occurs in the interior 
of S, it is a unique maximum and is given by 




( 2 . 12 ) 



for i = 1,2, 3, 4. where 
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0 < < 1 , 0 < ^; < 



for i = 2, 3, 4. 

In this case the MLE's for i = I, 2, 3, 4 are 



Qx 



fi 



(2.13) 



Xi 




«i 



(2.14) 



for i = 2, 3, 4. Finally we can estimate the reliability of the device as 



^ — QiQiQzQa- 



(2.15) 



Table 2 summarizes the cases for which the MLE's can be found using (2.12) - 
(2.15). When the maximum falls on the boundary of S, there is no explicit expression 
for the MLE of ( /?, , . . . , R^) . These are the cases with the exception of ( 1 1 2 1 ) 

A 

which have one failed item in the manufacturer test. This implies that R^ will be less 
than 1.0 . To find the MLE, we find the maximum of ( /?, , . . , R^) on each of the 
boundaries, compute / for each of these and let the MLE be the one with the largest 
value of / . It is clear that in most cases several of the boundaries can be eliminated 
from consideration, simplifying computation considerably. 
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Table 2. EASY AND HARD CASES 



EASY CASES 


HARD CASES 


( 


0 


0 


0 


0 ) 


( 1 


0 


0 


0 


) 


( 


0 


0 


0 


1 ) 


( 1 


0 


0 


1 


) 


( 


0 


0 


1 


0 ) 


( 1 


0 


1 


0 


) 


( 


0 


1 


0 


0 ) 


( 1 


1 


0 


0 


) 


( 


0 


0 


I 


I ) 


( I 


0 


I 


1 


) 


( 


0 


1 


0 


1 ) 


( 1 


I 


0 


1 


) 


( 


0 


1 


1 


0 ) 


( 1 


1 


I 


0 


) 


( 


0 


0 


2 


0 ) 


( I 


0 


2 


0 


) 


( 


0 


1 


1 


1 ) 


( 1 


I 


I 


1 


) 


( 


0 


0 


2 


» ) 


( 1 


0 


2 


I 


) 


( 


0 


I 


2 


0 ) 


( 1 


I 


2 


0 


) 


( 


0 


1 


2 


1 ) 




( 


1 


I 


2 


1 ) 





D. EXAMPLES 
1. Easy Case 

In this example our failure vector is, 



( 0 1 2 1 ). 



13 



According to this failure vector, we can write our likelihood equation by using equation 
(2.6) to get 

L = Rf{ 1 - R, )°Rl\ \ - R2? 1 - ^4 )' (2.16) 

0 < R^ < 1 , 0 <Ri <Ri , 

for i = 2, 3, 4. From the likelihood above; 






R2 



^3 



Ra 



^ " f- = '0000 

^ = W - 0-«oo 

= 0.9375 

«3 32 

•^4 19 

^ ^ = 0.9500. 



(2.17) 



A A 

Results imply that all R's for i = 2, 3, 4 are between 0.0 and R^, which means that 
constraints are met. Now we can estimate Q-s using equations (2.13) and (2.14). 



Qi = 


1.0000 






Qi = 


0.9500 

1 


= 0.9500 




^3 = 


0.9375 

1 


= 0.9375 


(2.18) 


Qa = 


0.9500 

1 


= 0.9500. 





And finally, reliability of the device can be estimated by using equation (2.15) 
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R = ( 1.0000) (0.9500) (0.9375) (0.9500) = 0.84609375. 



(2.19) 



2. Hard Case 

In this example our failure vector is 



(10 2 1 ) 



According to this failure vector, we can write our likelihood equation by using equation 
(2.6) as 



Z. = ( 1 - Ri )’ rI° ( 1 - /?2 )° ( 1 - /?3 )^ /?4’ ( 1 - /?4 )' (2.20) 

0 < /?, < 1 , 0 <Rj <Ri 



for i = 2, 3, 4. From the likelihood above; 



^1 19 

-»r = 

- -i- - i-ooeo 

— = — = 0.9375 
"3 32 



ii 

«4 



= 0.9500. 



( 2 . 21 ) 



As you see from the above thus, the MLE does not lie in the interior of S. 

"2 "I 

We begin by considering the boundary 
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/?] — /?2 — /?3 ~ ^ • 



The likelihood equation on this boundary is, 



L = ( 1 - 

0 < /?, < 1 . 



From the likelihood equation, 



A A A A CC 

/?, = = /?3 = = 0 . 9565 . 



Then value of the likelihood with these estimated R's is 



L = ( 



88 

92 



f ( = 0-714 X 



10 



-7 



For the boundary 



/?] — /?2 — ^3 > ^ ^ -^1 » 



the likelihood equation is. 



L = 1 - /?, 1 - /?4 )' 



0 </?, < 1 , 0 < R 4 < R^ . 



From the likelihood equation, 



A A A A 1 Q 

/?, = = /?3 = ^ = 0.9583 R ^= -^ = 0 . 9500 , 



( 2 . 22 ) 



( 2 . 23 ) 



( 2 . 24 ) 
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and the value of the likelihood with these estimated R/s is 



i = ( -ff- ^ ^ ^ = 0.723 X I0-’ 



For the boundary 



— R2 — R4 , R3 ^ R\ 



the likelihood equation is, 



L = 



I 



R, f ( 1 



- ^3 f 



0 < /?, < 1 0 < /?3 < /?, . 



From the likelihood equation. 



A A A <0 

/?, = = /?^ = ^ = 0.9666 



30 

32 



/?3 = -^ = 0.9375 . 



Then value of the likelihood with these estimated R's is 



^ )“< W )'< IT )”( ^ ">■’ 



For the boundary 



^1 ~ ^3 ~ ^4 » ^2 ^ ^1 » 



the likelihood equation is. 



(2.25) 



(2.26) 



(2.27) 
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(2.28) 



L = \ - Ry f rI° 

0 < Ry < I 0 < R 2 < Ry . 

Clearly L is maximized on the boundary of the constraints (2.28), i.e. on another of the 
boundaries of S thus we can eleminate this case from consideration. For the boundary 

Ry = R 2 , R 3 ^ R\ ) R/l ^ R\ t 

the likelihood equation is, 

L = /?r ( 1 - )' ^ 3 ° ( 1 - /ej )^ ( 1 - ^ )' (2.29) 

0 < Ri < 1 , 0 < ^3 <Ry , 0 < R4 <Ry. 



From the hkelihood equation. 



A A '70 ^ -^19 

i?, = /?2 = ^ = 0.9750 T?3 = ^ = 0.9375 ^ ^ = 0.9500 



Then value of the likelihood with these estimated R's is 



^ ^ ^ ^ ^ ^ 10 -’< 2 . 30 ) 



Quick inspection reveals that the remaining boundaries can be eliminated from consid- 
eration. 

After boundary analysis, we can see that equation (2.30) gives us the maximum 
likelihood value. Thus the MLE's for /?, , R 2 , R^, R^ are 

= 0.9750 
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y?2 = 0.9750 
^3 = 0.9375 



4 = 0.9500 . 



Now we can estimate Q,'s with equations (2.13) and (2.14) 



2, = 0.9750 



(23 

24 



0.9750 

0.9750 

0.9375 

0.9750 

0.9500 

0.9750 



1.0000 
0.9615 
0.9744 . 



(2.31) 



And finally reliability of the device can be estimated by using equation (2.15) 



R = ( 0.9750 )( 1.0000) (0.9615) (0.9744) = 0.91346100. (2.32) 



We computed MLE's with a FORTRAN program MLEA given in Appendix 
A. After 24 replications of this program, we can estimate R, 's for i = 1, 2, 3, 4 for each 
of the 24 outcomes of the sampling plan that lead to lot acceptance. These are given in 

A 

Table 3, R is given in Table 4. 

We compute lower confidence bounds for each possible case by bootstrapping. 
The bootstrap can be used to produce approximate confidence intervals in an automatic 
way. There are several ways to set approximate confidence intervals with bootstrapping. 
These are the percentile method, the standard method, the bias-corrected percentile 
method and the nonparametric method [Ref 3 ; pp. 67-70). Let d be an unknown pa- 

A A * 

rameter with estimator 6. To bootstrap, samples are generated using 6 in place of the 
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unknown parameter Q. For each of these samples an estimator O' is computed. We 

A ^ ... . . 

define (7 ( s ) to be the parametric bootstrap cumulative distribution function of 9', i.e. 

A ^ ^ 

(7 ( s ) is the empirical distribution of the 9‘'s. All methods mentioned above use 

A 

percentiles of G to define the confidence interval. They differ in which percentiles are 
used. The percentile method was used in our calculations. 

If we use the notation 6 {a) for the level a. endpoint of one sided lower confi- 
dence interval for 9 , then 



/>[0(a)<0] = a. (2.33) 

An estimate of 9 ( a ) from the bootstrap cumulative distribution is given by 



0(a)=G"‘(a). (2.34) 

To get LCB's for R, for each of the 24 realizations of the sampling plan , we 
generate bootstrap samples of random failure vectors, in which failures come from in- 

A 

dependent binomial distribution with parameters (n, , RJ for i = 1, 2, 3, 4. This is done 
using the FORTRAN program RANVEC [Ref. 4) in Appendix B. We generate 5000 

A 

failure vectors for each case. And then we estimate R 's from each of the 5000 failure 
vectors by the means of the program MLEA. The next step is to compute the order 

A A A 

statistics of R 's from R, to /?jooo. We get the parametric bootstrap cumulative distrib- 

A A 

ution function ( i.e. the emprical distribution of /?, , . . . , Rf„ao ) with this computation. 
This is done using the FORTRAN program SORT in Appendix C. 

Finally we can get the 95 % lower confidence bound using equation (2.34) from 
this routine. Reliabilities and 95 % lower confidence bounds are listed in Table 4. Re- 

A 

suits in Table 4 are given in descending order and R's and LCB's are not ordered as we 
expected them to be. For example, the failure vector that has the maximum number of 
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A 

failure in each test is in the middle of the table with respect to R. It has also a bigger 
95 % LCB than the failure vector that has a total of 3 failed items in each of the joint 
tests. These results are counter-intuitive because we expect that more failures indicate 
a lower overall reliability. This is reasonable because it is likely that a device that is 
poorly constructed is more likely to fail in any of the four environments. From the re- 
sults in Table 4, it is clear that an attempt to model dependence must be made in order 
to get believable estimates of R. 
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Table 3. TEST PROBABILITIES ( A ) 





A 


A 


A 




CASE 






R, 


k 


( 0000 ) 


1.0000000 


1.0000000 


1.0000000 


1.0000000 


( 0001 ) 


1.0000000 


1.0000000 


1.0000000 


0.9500000 


( 0010 ) 


1.0000000 


1.0000000 


0.9687500 


0.1.00000 


( 0100 ) 


1.0000000 


0.9500000 


1.0000000 


1.0000000 


( 1000 ) 


0.9891304 


0.9891304 


0.9891304 


0.9891304 


(001 1 ) 


1.0000000 


1.0000000 


0.9687500 


0.9500000 


( 0101 ) 


1.0000000 


0.9500000 


1.0000000 


0.9500000 


( 0110 ) 


1.0000000 


0.9500000 


1.9687500 


1.0000000 


( 1001 ) 


0.9861111 


0.9861111 


0.9861111 


0.9500000 


( 1010 ) 


0.9833333 


0.9833333 


0.9687500 


0.9833333 


( 1100 ) 


0.9861111 


0.9500000 


0.9861111 


0.9861111 


( 0020 ) 


1.0000000 


1.0000000 


0.9375000 


1.0000000 


( 0111 ) 


1.0000000 


0.9500000 


0.9687500 


0.9500000 


( 1011 ) 


0.9750000 


0.9750000 


0.9687500 


0.9500000 


( 1101 ) 


0.9807692 


0.9500000 


0.9807692 


0.9500000 


( 1110 ) 


0.9750000 


0.9500000 


0.9687500 


0.9750000 


( 0021 ) 


1.0000000 


1.0000000 


0.9375000 


0.9500000 


( 0120 ) 


1.0000000 


0.9500000 


0.9375000 


1.0000000 


( 1020 ) 


0.9833333 


0.9833333 


0.9375000 


0.9833333 


( 1111 ) 


0.9615384 


0.9500000 


0.9615384 


0.9500000 


( 0121 ) 


1.0000000 


0.9500000 


0.9375000 


0.9500000 


( 1021 ) 


0.9750000 


0.9750000 


0.9375000 


0.9500000 


( 1120 ) 


0.9750000 


0.9500000 


0.9375000 


0.9750000 


( 1121 ) 


0.9500000 


0.9500000 


0.9375000 


0.9500000 
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Table 4. RELIABILITIES AND 95 % LOWER CONFIDENCE BOUNDS 



(A) 



FAILURE 

VECTOR 


A 

R 


FAILURE 

VECTOR 


95 % LCB 


(0000) 


1.0000000 


(0000) 


1.0000000 


(1000) 


0.9891304 


(0010) 


0.9062500 


(0010) 


0.9687500 


(1000) 


0.9000000 


(1010) 


0.9687499 


(0020) 


0.8750000 


(0001) 


0.9500000 


(1010) 


0.8550000 


(0100) 


0.9500000 


(0001) 


0.8500000 


(1001) 


0.9499999 


(0100) 


0.8500000 


(1100) 


0.9499999 


(1001) 


0.8282812 


(1011) 


0.9439102 


(1100) 


0.8258822 


(1110) 


0.9439102 


(01 10) 


0.8234775 


(1111) 


0.9385999 


(0011) 


0.8234375 


(0020) 


0.9375000 


(1020) 


0.8125000 


(1020) 


0.9374999 


(1011) 


0.8015624 


(1121) 


0.937499 


(0101) 


0.8000000 


(0011) 


0.9203125 


(1110) 


0.8000000 


(0110) 


0.9203125 


(0120) 


0.7875000 


(1101) 


0.9201959 


(0021) 


0.7749999 


(1021) 


0.9134614 


(1101) 


0.7749018 


(1120) 


0.9134614 


(1111) 


0.7649999 


(0101) 


0.9025000 


(1021) 


0.7647058 


(0021) 


0.8906249 


(1120) 


0.7614843 


(0120) 


0.8906249 


(1121) 


0.7505192 


(0111) 


0.8742968 


(0111) 


0.7499999 


(0121) 


0.8460937 


(0121) 


0.7171874 
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III. LOG LINEAR MODEL WITH DEPENDENCE ASSUMPTION 



A. BACKGROUND 

Modeling the outcomes of the various environmental tests as independent is clearly 
inappropriate. However, the nature of the acceptance sampling plan makes it impossible 
to estimate the reliability of the device ( i.e. the probability that it would activate after 
exposure to all four environments ) without some assumptions about the dependence 
between outcomes of various tests. Thus, our second approach is to model the 
pyrotechnic device reliability with a log linear model. 

The results of our test series create a(2x2x2x2) contingency table. Let p 
represent passing a test and f represent failing a test. 



Table 5. CONTINGENCY TABLE STRUCTURE FOR TEST SERIES 





Passed 

Manufacturer Test 


Failed 

Manufacturer Test 


Passed 

Tem.&Hum. 

Test 


Failed 

Tem.«S:Huni. 

Test 


Passed 

Tem.&Hum. 

Test 


Failed 

Tem.&Hum. 

Test 


Passed 

Vibration 

Test 


Passed 

Altitude 

Test 


^pppp 


^Pfpp 


^fPPP 


^^ffpp 


Failed 

Altitude 

Test 


^PPPf 


^^pfpf 


^fppf 


^Jfpf 


Failed 

Vibration 

Test 


Passed 

Altitude 

Test 


^ppfp 


^^pffp 


^^fpfp 


^^Mp 


Failed 

Altitude 

Test 


^ppff 


^pJJT 


^fpJT 
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The results of our tests can be thought of as censored data from a hypothetical 
(2X2X2X2) contingency table (See Table 1) . The frequency in each cell of this table 
is 



{ij, k, 1) e {p, ff , 



the number of devices out of 92 which would have result i in environment 1 (manufac- 
turer), result j in environment 2 (temperature and humidity alone), result k in environ- 
ment 3 (vibration alone), result I in environment 4 (altitude alone). 

This is a hypothetical table, because if a device was exposed to all four environments 
and then failed to activate, there would be no way to discern which combination of the 
four environments caused failure. The data from the acceptance sampling plan can be 
thought of as censored data from such a (2X2X2X2) contingency table. As an 
example, a device that is given just the manufacturer test belongs in one of the cells 
( p, j, k, / ) where ( j, k, / ) e {p, /}^ because it is not clear what would have happened 
to it had it been exposed to the other three environments. 

Using a log linear model under independence, the expected value of each cell fre- 
quency is 






(3.1) 



where i, j, k, / = p, f and 






i 

P 




(3.2) 



for i = 1,2, 3, 4. Here the parameter p represents the overall effect and the parameters 
2' represent the effect of passing environment i. It is simple to show that the log linear 
model (3.1) is equivalent to the independence of the environments [Ref 5; pp. 25-46]. 
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On a logarithmic scale, the independence relation is equivalent to the additive re- 



lationship. As an example, 



log ^pppp = + ^p + 



(3.3) 



Log linear models, which take into account dependence include extra interaction terms. 
Because of the extreme amount of censoring, we will only consider models with two way 
interaction terms of the type. 



_ ll . l'2 . ,13 , ,14 , ,23 , ,24 ,34 



Hjkl 



(3.4) 



where i, j, k, / = p, f, 



4 = - 4 



(3.5) 



for i = 1, 2, 3, 4 and 



4 = 4 = -4f= -4> 



(3.6) 



for i, j, = 1,2, 3, 4. As an example; 



E[ M 1 = 

^ L ^»pppp J c 



24 . *34 



(3.7) 



The censored data does not tell us much about partial association between pairs of 
tests. For this reason, we assume that partial associations between any pair of tests are 
all the same. Let 6 be the partial association between tests. Then we can reformulate 
equation (3.7) as follows: 



E C Mpppp ] 



+ X^ + X* + 6 0 



(3.8) 
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The MLE's for the expected number of devices in each cell ( or equivalently the 
MLE's for the parameters ) can not be found explicitly. The EM Algorithm ( Expecta- 
tion - Maximization ) will be used to approximate the MLE's of the expected number 
of devices in each cell. 

B. EXPECTATION MAXIMIZATION ( EM ) ALGORITHM 

A general method of maximum likelihood estimation from incomplete data is the 
E.M Algorithm. The expectation maximization algorithm is an iterative procedure where 
each stage consists of: 

■ an expectation step ( E ) followed by 

■ a maximization step ( M ). 

This algorithm is generally used to compute maximum likelihood estimators in in- 
complete data problems. In the application of the EM algorithm we: 

• replace missing values by their estimated expected values given the incomplete data 

• estimate parameters 

• reestimate the missing values assuming the new parameter estimates are correct 

• reestimate parameters 

and so forth, iterating until convergence (Ref. 6: pp. 127-141]. This iterative algorithm 
works as follows in pyrotechnic device problem. 

1. Initialization 

The algorithm requires initial guesses for the parameters of the log linear model 
cell frequency in contingency table. It uses initial guesses and calculates initial cell 
probabilities. We will get our initial guesses by first fitting a log linear model under in- 
dependence. 
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2. Iterations 



■ Expectation Step 

EM Algorithm estimates expected cell frequencies given the data by using the 
most current estimates of the parameters. It compares these estimated expected cell 
frequencies with the previous estimates. If the differences between these two estimates 
are small enough, then the EM Algorithm has converged. We then accept the final es- 
timates of cell frequencies as MLE, so we can easily calculate cell probabilities. 

■ Maximization Step 

If the algorithm has not converged, then it starts to estimate new parameters for 
the log linear model using expected cell frequencies from the E step as if there were ac- 
tual data available. Estimation is done by maximizing the likeUhood using an iterative 
Newton- Raphson method. The estimated parameters from this step are than used in the 
next E step of the EM algorithm. 

C. CALCULATIONS 

We apply the EM algorithm for each realization of the failure vector ( FOM, 
FOTH, FOV, FOA ). During the maximization step of the EM Algorithm, we use the 
Newton-Raphson procedure which is described by SAS, (Ref 7: pp. 190-212). We first 
describe the M step using as an example the following failure vector. 

(112 0 ) 

■ Maximization Step 

We start with initial guesses (IG) for A/^*/s. These are the conditional expected 
values provided by the previous E step. For example, initial guesses for the failure vector 
above, are shown below. 
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IG 



^hppp 

‘^fpppf 



M, 



ppfp 



M 



ppff 



M 



pfpp 



M 



pfpf 



M 



pffp 



M 



pfff 



^^fppp 

^'Ifppf 

^^fpfp 

^^fpff 



^tffpp 

^^fffP 

^ffff 



IG 



83.08 

0.1 

3.36 

0.1 

2.21 

0.1 

0.1 

0.1 

2.15 

0.1 

0.1 

0.1 

0.1 

0.1 

0.1 

0.1 



( 3 . 9 ) 



From these, we can compute the proportion of observations which fall into each cell as; 



P 



M, 



pppp 



92 



M, 



pppf 



92 



M, 



ppfp 



92 



0.9030 

0.0011 

0.0365 



( 3 . 10 ) 



92 



0.0011 



A 

We will use P to get MLE's for the parameters of the log linear model. Our log linear 
model which includes two way interactions can be written as 
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(3.11) 



log Ppppp 
log Ppppf 
log Pppfp 

log Pffff 



’ll 1 1 16' 

11 1 1-10 


y - p 


* 1 ' 
1 
1 


1 -1 -1 -1 -1 6 




T 

t 



where = {n, A), A* ) is the parameter vector and p is normalizing constant 

required by the restriction that the probabilities sum to 1 and where is the proba- 
bility corresponding to the {ij klY cell in the hypothetical contingency table. We can 
derive from equation (3.1 1) 



log Ppppp 
log Ppppf 
log Pppfp 




' 1 1 1 16’ 

1 1 1-10 


^ - d 


’ 1 ' 
1 
1 


\OgPffff 




-1 -1 -1 -1 6 




1 



where = { A^, A^, /IJ, A* ) is the parameter vector and <5 is normalizing constant re- 
quired by the restriction that the probabilities sum to 1. Then, to use the SAS proce- 
dure, we rewrite log log . . . ,logP^ as the 15 logits, F, = ( log P^J log Pff„), 
. . . , = ( log Pfff^f log Pffff) so that 



F{P) 



' Fx ' 




’ 1 0 0 0 -1 ■ 


F2 





0100-1 


Fx5 

<- - 




0000-1 



log Ppppp 
log Ppppf 
log Pppfp 

^OgPffff 



(3.13) 
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where P = ( , Pffff). Using equations above the following result is ob- 

tained 



F{P) 



‘ Px ' 




' 2 2 2 2 o' 


Pi 





2220-6 


Fxs 




0002-6 



( 3 . 14 ) 



The design matrix (X) of this problem, from the equation (3.14) is as follows: 



X 



2 2 2 2 0 
2220-6 
2202-6 
2200-8 
2022-6 
2020-8 
2002-8 
2000-6 
0222-6 
0220-8 
0202-8 
0200-6 
0022-8 
0020-6 
0002-6 



(3.15) 



In the application of Newton-Raphson Method, we use the variance and covariance 

A 

matrix S of F( /* ). It's inverse is given by 



s-'(P) 



P 

^pppp 


P^ 

‘ pppp 


— P P 

PPPP PPPf , 


— Ppppp ^ Pfffp 




_ P 

^PPPP 


^ ^pppf 


P — P^ 

PPPf PPPf 


• • ~ ^PPPf ^ ^fffP 


.(3.16) 






. . 


- P 
^PPPP 


^ ^fffp 


~ ^PPPf ^ ^fffp 


1 
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The first estimates of parameters for the Newton- Raphson procedure are calculated 



as follow's: 






( 3 . 17 ) 



where 5~‘ and F are estimated using the proportions in the vector P. We estimate the 
reduced logit response functions using the equation 



A 



Fo = Xb,. 



( 3 . 18 ) 



From Fo we can compute the updated estimates of P„,f, •••, Pfffj Let 

n = ( Hi ( 1 ), 03(1), . . . , n,s ( 1 ) ) be the vector, which contains these estimates of 
the cell probabilities. The value of the log likelihood evaluated at fl is 



where x, is the nurriber of items in cell i with respect to probability in Fl. 

We estimate parameters iteratively until the difference between last estimate and 
previous is small enough. At each iteration we update the inverse of variance and 
covariance matrix W’ith probabilities of fl from the previous iteration and we do fol- 
lowing matrix computations. 



Let b, be the next estimate of parameters /? in the i th iteration. Then F follows: 



16 




( 3 . 19 ) 



c = (n)x 

G = X^[92.0 X (F - n)] 



( 3 . 20 ) 



bi = bi_^ - SC~^G 



( 3 . 21 ) 
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where 5 < 1 is a constant supplied by the user. 

We get first cell expectations from initialization. The failure vector, which is defined 
at page 5, has initial guesses for cells which are shown in equation (3.9). After the first 
application of Xewlon -Raphson Algorithm, we have following initial cell expectations 
(CE ). 



82.4278259 

0.7278117 

3.1047230 

0.1208856 

2.0790358 

0.0809494 

0.3453168 

0.0592893 

2.0268116 

0.0789160 

0.3366428 

0.0578000 

0.2254283 

0.0387050 

0.1651089 

0.1250073 



(3.22) 



■ Expectation Step 

The conditional expected frequency for each cell expectations are calculated using 
the estimated parameters of the log linear model from the previous M step. Some of the 
cell expectations formulated below 



E iMpppp I {FOM, FOTH, FOV, FOA) ] = (20 - FOM) xP{pppp\p...) 



+ (20 - FOTH) X P ip p p p \ p p . 
+ (32 - FOV)xP{pppp\p.p.) 
+ (20 - FOA) xP{pppp\p..p) 
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E\_Mpjpf\{FOM,FOTH,FOV,FOA):\ = (20 - FOM)x P {pfpf\ p . . .) 

+ {FOTH)xP{pfpf\pf..) 

+ (32 - FOV)xP(pfpf\p.f.) 
+ (FOA) X P{pfpf\p..J) 



where " . " represents either a pass or fail. For example, 

F{pppp\p---)is the probability that device passes all environments given that it 
passed manufacturer test. 

P {p p p p \ p p . .) is the probability that a device passes all environments given that 
it passed manufacturer test and temperature-humidity environment. 

The conditional probabilities are computed from the estimated cell probabilities 
from the previous M step. As an example, we have following failure vector, 

(112 0 ) 



Then F C | ( 1, 1, 2, 0) ] is calculated as follows: 



P{PPPP\P‘^-) 



82.4278259 

88.9458375 



0.9267193 



n/ . ^ 82.4278259 ^ 

P(pppp\pp..) = = 0.9542328 



n, . . 82.4278259 

P(pppp\p.p.) - - 0.9661516 



n/ 1 ^ 82.4278259 ^^.^..,00 

P(PPPP\P-P)- 87,,5^,o15 -0«71388. 



where the numbers come from (j.22) the estimated cell expected values after from the 
previous M step. Thus, 
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ElMpppp\{\, 1, 2, 0)] = 83.4654163 



The remaining 16 conditional expectations are calculated similarly. Then we com- 
pare these expectations with previous expectations. If the difference between compared 
expectations are small enough then EM Algorithm is assumed to have converged. 

When EM Algorithm converges, we can estimate reliability of device ( R ) after 
exposure all four environments as 



R = 



^ [ ^pppp ] 



92.0 



(3.25) 



For example, for the failure vector above, after 27 iterations, we obtain 

E\_Mpppp'] = 84.412714, 



and 



R = 



84.412714 

92.0 



0.9175259. 



D. INITIAL GUESS PROBLEM 

The EM Algorithm uses an initial guess vector as in equation (3.9). We use the 
independence assumption as described in equations (3.1), (3.2) and (3.3). The following 
procedures are used during the application of EM Algorithm to the log linear model with 
the independence assumption. 

We have still an initial guesses problem for parameters in the log linear model with 
the independence assumption. For this reason, we do some calculations to find initial 
guesses for the parameters. Moreover, one can choose random initial guesses for the 
parameters too. 
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Assume that we have the following failure vector, 



(10 10 ). 



In this procedure we divide the number of devices which pass the manufacturer test 
by the number of cells which include items that have passed the manufacturer test. In 
this example, we can calculate as follows: 















n 3 = 2.375 items for manufacturer test 

o.O 

t 3 = 2.500 items for temperature and humidity test 

8.0 

C 4 ^ 3 = 0.125 item for vibration test 
8.0 

C 3 = 0.000 item for altitude test 
8.0 



fmally the expected number of device for this cell is. 



^ppff ~ • 



After computing for each cell using the same procedure above, we can group these 
cells to estimate , A/,^ as follows: 

* — ^^ppOO 

* ^^p/kl ~ ^p/OO 

* Yi^fpkl ~ ^fyOO 

* ^^//kl ~ 

where k, / = p, f. 
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In this example. 



Table 6. GRUOPED DATA FOR INITAL GUESS CALCULATION 



TERMS 


ESTIMATED VALUES 


^fppOO 


32.500 




22.500 




23.500 




13.500 



By combining the parameters and XI , we have four equations and four unknown 
parameters. We can easily solve for the parameters from the initial guesses. 



Table 7. ESTIMATED PARAMETERS FOR INITIAL GUESS 



TERMS 


ESTIMATED VALUES 


/i .+ K 


-1.97 


K 


1.08 


K 


1.08 


K 


0.56 



Finally, we can split the sum of the n and X^ into two parts. One possibility is to 
divide the sum by two assigning half to ^ and half to Xj,. Results using a random initial 
guess and the above procedure for an initial guess are very close. Calculation of initial 
guess parameters are done by the FORTRAN program INITIAL in Appendix D and 
the FORTRAN program PARAM in Appendix E. 
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E. RESULTS 



Reliabilities of the pyrotechnic device are calculated by a FORTRAN program 
LLMDEP in Appendix F. They are given in Table 8. Results from the log linear model 
design with dependence assumption are similar to the first model. We can easily see 
from the previous table that the ordering of the estimates of R is counter-intuitive. We 

A A 

expect ( 1 1 2 1 ) to yield the smallest R rather than an R = 0.934 which is larger than 

A 

R for about half of the 24 realizations of the failure vector. 

It is clear from these results that this log linear model is inappropriate for modeling 
the outcomes of the sampling inspection plan. More realistic models would include 
three-way and four- way interaction terms. However, due to the extreme censoring in the 
data, we can not estimate R for these models. Moreover, reliabilities of several cases for 
even this model were not calculated because of computational limitations. Taking the 
account of dependence with log linear models is clearly not a reasonable approach for 
estimating reliabilities from this data. 
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Table 8. RELIABILITIES WITH LOGLINEAR MODEL 



FAILURE VECTOR 


R 


(0000) 


1.0000000 


(1000) 


0.9887247 


(1010) 


0.9693505 


(1001) 


0.9518477 


(1100) 


0.9503552 


(1110) 


0.9457148 


(1011) 


0.9456805 


(1111) 


0.9399773 


(1020) 


0.9374950 


(1121) 


0.9336501 


(0011) 


0.9255090 


(0110) 


0.9255087 


(1101) 


0.9246104 


(1021) 


0.9175486 


(1120) 


0.9175295 


(0101) 


0.9089977 


(0021) 


0.8978485 


(0120) 


0.8978109 


(0111) 


0.8825021 


(0121) 


0.8584540 
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IV. WORST CASE SCENARIO 



A. ASSOCIATION ANALYSIS 

In the inspection sampling plan, failures are not assigned to an individual failure 
mechanism except for the manufacturer test. For example, from the following failure 
vector. 



(1 12 1 ) 

We know that there is at least one failure due to the manufacturing test. There is one 
failure from the joint temperature and humidity and manufacturer test, but we do not 
know which failure mechanism generated this failure. Items can fail due to a manufac- 
turing related failure mechanism or other failure mechanism or both. The worst case is 
to assume that the cause of failure is due to all of the failure mechanisms that it was 
exposed to. 

Let MT represents manufacturer test, TAHT represents temperature and humidity 
test, VT represents vibration test, AT represents altitude test. 

Because all 92 items are exposed to manufacturer test, the worst case of the failure 
analysis of the failure vector above is shown by the following table. 



Table 9. FAILURE ANALYSIS 





MT 


TAHT 


VT 


AT 


# OF TESTED ITEMS 


92 


20 


32 


20 


# OF SUCCESFUL 
ITEMS 


87 


19 


30 


19 



A 



This " worst case " scenario should give us lower bounds for the estimates R of R. 



If vve calculate the reliability of scries system assuming failure mechanisms in test 
are independent when in fact they are associated but not independent, then we underes- 
timate system reliability. It is reasonable to expect that there is positive dependence 
between four tests. If one item fails in manufacturer test, it is more likely to fail in other 
tests. One way of modelling positive dependance between tests is to assume that tests 
are positively quadrant dependent. 

Given random variables T^ T^. They are said to be positively quadrant de- 

pendent (PQD) [Ref. 8: p. 33], if 

n 

P(r, < q, ..., < r„) > ^ h) (4.1) 

i=i 

for all ( r,, r^, . . . , r„ ) e R". An equivalent formulation of positive quadrant depend- 
ence is r,, Ti, . . . , V' are PQD iff 

n 

P{T,> t,, t„)^\\p(T,> ti) (4.2) 

j=i 

for all ( r„ . . . , r, ) e /?". The proof that (4.1) and (4.2) are equivalent as given in 
[Ref. 8: pp. 32-33]. 

We may take account positive dependence between tests with the PQD assumption. 
Let 



Ti = 



1 , if an item passes environment i 

0 , if an item fails envirorunent i 



(4.3) 



for i = 1,2, 3, 4. Assume that T, , T^, Tj , T^, are positively quadrant independent, i.e. 
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(4.4) 



<4 

P ( 7', < /] , T 2 < i 2 , T 2 < I 2 , < t^) > {Tj < tj ) 



/=i 



or equivalently that 



P(r, > t,, T 2 > l 2 T 2 > I 3 , > u) > > /,). (4.5) 

1=1 



Let R be the probability that an item activates after exposure to all 4 environments. 



Then 



/? = P( r, = 1 , Tj = 1 , T3 = 1 , T, = 1 ) . (4.6) 

Using equations (4.4) and (4.5) 

R> P( r, = 1)P(T2= 1)P{T2 = l)P{T^ = 1 ) . (4.7) 

Using the notation from the previous section 

(4.8) 

where is the probability that an item passes environment i. With the censored data, 
we can't estimate Q^ , Q^, and Q, without building a more structured model for 
r, , T 2 , Ty and T^. One alternative is get a lower bound for estimates of 
Qit Qit Qa' This is worst case scenario . 

B. CALCULATIONS WITH EXAMPLE 
Let 

■ R be the MLE for the reliability of the device. 

■ Qi be the MLE for the probability that an item passes environment i. 



42 



■ /? be a lower bound of R for reliability of the device. 

A 

■ be a lower bound for the probability that an item passes environment i. 
Then according to equation (4.8) 



R > 




(4.9) 



Further we will construct estimates Q, of Q, such that 



4 ^ 4 ( 4 . 10 ) 



and then define 



Thus 




(4.11) 



R < L (4.12) 
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From Chapter II, the likelihood of observing JT, = x,, A’j = Xj, A'j = Xj, is 



» *^ 2 ’ *^ 3 * *^4 1 ^ 2 * ^ 3 * ^4 ) 



4 







p(7, = 1)^' P(r, = o)”'"'^' 



X = 1.72=1)"^ 

X (1 -P(7i = 1,72= (4.13) 



X 7(7, = 1,73=1)"^ 

X (1 -7(7, = 1 , 73 = 1) 

X 1 ^""^ 7 (7, = 1 , 74 = 1 )^ 

X (1 -7(7, = 1,7»= 1))"^-^ 



where x, is the number of devices out of «, that activate after test i. If we know why the 
device failed for the tests i = 2, 3, 4 which include manufacturer test along wdth exposure 
to environment i then our likelihood could be written as 
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L(jr,, jTj, X3, x^\R^, R2, /?3, ^4) = 




p(r, = i)^> p(r, =o)'’> 



jti 



X 



/J2 



^•^ 2 » -^ 21 * '^ 22 * '^ 23 / 



X p (r, = 1 , T2 = P (r, = 1 , 72 = o)^^‘ 

X 7 ( 7 ,= 0 , 72 = 7 ( 7 , = 0 , 72 = 



X 



■' ) 

-■^3> -^31- •^32* •^33/ 



(4.14) 



X 7 ( 7 , = 1 , 73 = 7 ( 7 , = 1 , 73 = 0 )^^‘ 

X 7 ( 7 ,= 0 , 73 = 7 ( 7 , = 0 , 73 = 0 )''” 



X 



-^41’ *^42> *^43> 



X 7 ( 7 , = 1 , 74 = 1 )"^ 7 ( 7 , = 1 , 7 , = 0 )^‘ 
X 7 ( 7 ,= 0 , 74 = 1 )^^ 7 ( 7 , = 0 , 74 = 0 )^^ 



where jr,, is the number given test i which failed due to environment i but passed manu- 
facturing, is the number given test i that passed environment i but failed manufac- 
turing and x,j is the number given test i that failed both manufacturing and environment 
i for i = 2, 3, 4. Note that n, — x, = -f jc^ + x,j for i = 2, 3, 4. 

From this likelihood, the MLE's of <2; = P(P/= V are given by 



C, = 



Xj + (X2 + X21) + (jT3-bjr3,) + (jr4-bX4,) 

rt, + «2 + ^*3 + ”4 



( 4 . 15 ) 
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(4.16) 



Qi = 



•^2 ■^22 
n-^ 



ft = 



•>^3 + ^32 
«3 



ft = 



+ X42 



However since we do not know the x,/s, we see that 






jr, + jcj + 0:3 + X4 

n 



and 



ft ^ 



ii 



for i = 2, 3, 4. Let 



e, 



X| + ^2 + X3 + X4 

n 



and 




for i = 2, 3, 4. Then Qis are lower bounds for the true MLE's Q„ and 

A A A A A 

^ = 2i ft ft ft 

It is clear that ^ is a lower bound for the MLE R. 



(4.17) 



(4.18) 



(4.19) 



(4.20) 



(4.21) 



(4.22) 



(4.23) 
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In Worst Case Scenario, we assume that the item, which fails in test i, fails because 
of both manufacturing related failure mechanism and other failure mechanism which is 
induced by the i th test environment. With this assumption, the reliability of of device 
is given as follows 



R = Q^Q2 <23 <24 (4.24) 



where 



Qi 

A 

Q2 

A 

a 

A 

04 



JT, + Jf 2 + JC 3 + Jf4 

92 



20 

32 

ii 

20 



Here is an example, in this example our failure vector is, 



( 0 1 2 1 ) 



(4.25) 



4 = — = 0.9456521 

4 = ^ = 0.9500000 

4 = ^ = 0.9375000 

4 = ^ = 0.9500000 



(4.26) 



And finally reliability of the device can be estimated by using equation (4.24) 



R = ( 0.9456521 ) ( 0.9500000 ) ( 0.9375000 ) ( 0.9500000 ) = 0.8093069 (4.27) 
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C. RESULTS 



We compute MLE's with a FORTRAN program MLEB in Appendix G. These 
are given Tables 10 and 1 1. 

To get an approximate lower confidence bounds for R using the worst case data, 
we bootstrap using the procedure described in the previous chapter. The FORTILAN 
program RANVEC in Appendix B is used to generate the 5000 bootstrap samples for 
each case. We then estimate R 's for 5000 failure vectors by the means of MLEB for 

A A A 

given case. The next step is to compute the order statistics of R 's from R^ to Rjooo- This 
is done by FORTRAN program SORT in Appendix C. Finally, we obtain the 95 % 
lower confidence bound from this routine. Reliabilities and 95 % lower confidence 
bounds are tabulated in the following pages. Results of reliabilities and 95 % lower 
confidence bounds are given in descending order. 

These estimates and LCB's for R are decreasing with the total number of the failures 
out of 92 items tested. This, at least, is consistent with how w’e believe pyrotechnic de- 
vices behave. It should be noted that these estimates are in fact conservative lower 
bounds for the time MLE's under the very weak assumption of PQD. How' conservative 
these estimates are cannot be determined without more extensive data that allows us to 
estimate the degree of dependence between tests. 
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Table 10. TEST PROBABILITIES ( B ) 





' ' ' A 


A 


A 


' ' A ' ' 


CASE 




Rz 


Rz 


R. 


( 0000 ) 


1.0000000 


1.0000000 


1.0000000 


1.0000000 


( 1000 ) 


0.9891340 


1.0000000 


1.0000000 


1.0000000 


( 0010 ) 


0.9891340 


1.0000000 


0.9687500 


0.1.00000 


( 1010 ) 


0.9782608 


1.0000000 


0.9687500 


1.0000000 


( 0001 ) 


0.9891304 


1.0000000 


1.0000000 


0.9500000 


( 0100 ) 


0.9891340 


0.9500000 


1.0000000 


1.0000000 


( 1001 ) 


0.9782608 


1.0000000 


1.0000000 


0.9500000 


( 1100 ) 


0.9782608 


0.9500000 


1.0000000 


1.0000000 


( 0020 ) 


0.9782608 


1.0000000 


0.9375000 


1.0000000 


( 1020 ) 


0.9673913 


1.0000000 


0.9375000 


1.0000000 


( 0011 ) 


0.9782608 


1.0000000 


0.9687500 


0.9500000 


( 0110 ) 


0.9782608 


0.9500000 


0.9687500 


1.0000000 


( 1011 ) 


0.9673913 


1.0000000 


0.9687500 


0.9500000 


( 1110 ) 


0.9673913 


0.9500000 


0.9687500 


1.0000000 


( 0101 ) 


0.9782608 


0.9500000 


1.0000000 


0.9500000 


( 1101 ) 


0.9673913 


0.9500000 


1.0000000 


0.9500000 


( 0021 ) 


0.9673913 


1.0000000 


0.9375000 


0.9500000 


( 0120 ) 


0.9673913 


0.9500000 


0.9375000 


1.0000000 


( 1021 ) 


0.9565217 


1.0000000 


0.9375000 


0.9500000 


( 1120 ) 


0.9565217 


0.9500000 


0.9375000 


1.0000000 


( 0111 ) 


0.9673913 


0.9500000 


0.9687500 


0.9500000 


( 1111 ) 


0.9565217 


0.9500000 


0.9687500 


0.9500000 


( 0121 ) 


0.9565217 


0.9500000 


0.9375000 


0.9500000 


( 1121 ) 


0.9456521 


0.9500000 


0.9375000 


0.9500000 



49 



Table 11. RELAIBILITIES AND 95 % LOWER CONFIDENCE 



BOUNDS ( B ) 



FAILURE 

VECTOR 


A 

R 


FAILURE 

VECTOR 


95 % LCB 


(0000) 


1.0000000 


(0000) 


1.0000000 


(1000) 


0.9891304 


(1000) 


0.9891304 


(0010) 


0.9582200 


(0010) 


0.8766983 


(1010) 


0.9476901 


(1010) 


0.8766983 


(0001) 


0.9396738 


(0001) 


0.8222825 


(0100) 


0.9396738 


(0001) 


0.8222825 


(1001) 


0.9293478 


(0100) 


0.8222825 


(1100) 


0.9293478 


(1 10 0) 


0.8222825 


(0020) 


0.9171195 


(0020) 


0.8179346 


(1020) 


0.9069293 


(1020) 


0.8084239 


(0011) 


0.9003056 


(0011) 


0.7712974 


(0110) 


0.9003056 


(0110) 


0.7712974 


(1011) 


0.8903022 


(1011) 


0.7712974 


(1110) 


0.8903022 


(0101) 


0.7712974 


(0101) 


0.8828803 


(0101) 


0.7548369 


(1101) 


0.8730705 


(1101) 


0.7548369 


(0120) 


0.8615828 


(0120) 


0.7218070 


(0021) 


0.8615828 


(0021) 


0.7200747 


(1021) 


0.8519021 


(1021) 


0.7200747 


(1120) 


0.8519021 


(1120) 


0.7200747 


(0111) 


0.8457870 


(0111) 


0.6929346 


(1111) 


0.8362838 


(1111) 


0.6927614 


(0121) 


0.8093069 


(0121) 


0.6508338 


(1121) 


0.8001103 


(1121) 


0.6505433 
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V. BONUS SYSTEM APPROACHES 



A. BACKGROUND 

Quality is described as " especially high degree of goodness or worth " [Ref. 9 
: p. 685). In industry', a quality product is one that fulfills customer expectations. 
There are two general aspects of quality: 

■ Quality of Design 

■ Quality of Conformance. 

All goods and services are produced in various grades or levels of quality. These 
variations in grades or levels of quality are intentional; therefore the appropriate tech- 
nical term is quality of design. The quality of conformance is how well the product 
conforms to the specifications and tolerances required by the design. Quality of con- 
formance is influenced by a number of the following factors: 

■ the choice of manufacturing process 

■ the training and supervision of the workforce 

■ type of quality assurance system ( process controls, tests, inspections, etc. ) used 

■ the extent to which these quality-assurance procedures are followed 

■ the motivation of workforce to achieve quality. 

Quality Control is the engineering and management activity by which we measure 
the quality characteristics of a product, comparing them with specifications or require- 
ments and taking appropriate remedial action whenever there is a difference between the 
actual performance and the standard [Ref 10 : pp. 1-3). 

As contracts are now written, contractors need only to satisfy the requirments of the 
sampling inspection plan for lot acceptance. Contractors have no incentive to improve 
the quality of the items they provide, although they are in a position to do so. As 
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mentioned before, the quality of conformance is influenced by the motivation. Therefore 
we can motivate manufacturers by giving a bonus for improved quality. To improve 
quality, The Naval Weapons Support Center has decided to implement a bonus system 
for pyrotechnic devices. 

B. BONUS PLANS 

In this chapter, we design a bonus system to improve the quality of pyrotechnic 
devices. A good bonus system encourages reliability growth, because firms tr>^ to reach 
a high quality to get a bonus. An effective bonus system must detect small differences 
among the offered lots. In the previous chapters, we estimated the reliability of a 
pyrotechnic device in several w’ays. We assumed independence in the first model, and 
we assumed dependence betw'een tests in the last two models. The estimated reliabilities 
in the first two models are close to each other, but they exhibit different structures in 
order. One can easily see that the failure vector of the worst case (112 1) has bigger 
reliability value than ten of the possible cases, from Table 5 and Table 15. In addition, 
the order of 95 Vo LCB's of cases does not match to the order of MLE's in independent 
models. But in the worst case scenario, which assumes dependence between tests, we 
get the same order for both 95 % LCB's and MLE's. Thus, we will use the worst case 
scenario model in further calculations for the bonus system. With the assumptions 
above, we can apply three sampling plans for giving bonuses to the manufacturers. 
These are; 

■ Single Sampling Bonus System 

■ Double Sampling Bonus System 

■ Multi-Sampling Bonus System. 

Manufacturers have to meet pyrotechnic device acceptance criteria first, before 
having a chance to get a bonus. 
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1. Single Sampling Bonus System 

In this bonus system, we will decide two things at the end of inspection. First 
we will decide whether the offered lot is acceptable or not, and then for lots which are 
accepted we will decide whether to give a bonus to the manufacturer. 

The Single sampling bonus system is designated by three numbers. These are, 

M, , LCBFB , LCB 



where 

■ rt, is the sample size for inspection ( here n^ = 92 ). 

■ LCBFB is the cut off value applied to the lower confidence bound for 
awarding a bonus. 

■ LCB means estimated lower confidence bound after inspection. (If LCB is 
greater than or equal to LCBFB, then a bonus is aw'arded.) 

The following algorithm shows us how' single sampling bonus system works for 
pyrotechnic devices. 



STEP # 1 
STEP # 2 
STEP it 3 
STEP it 4 
STEP it 5 

STEP it 6 
STEP it 7 



DETERMINE CUT OFF VALUE FOR AWARDING THE BONUS. 
TAKE A SAMPLE SIZE OF 92 FROM OFFERED LOT. 

APPLY MANUFACTURER AND THREE ENVIRONMENT TESTS. 
COMPARE RESULTS OF TESTS WITH ACCEPTANCE CRITERIA. 
ESTIMATE ITS LOWER CONFIDENCE BOUND ( LCB ), 

IF THEY MET ACCEPTANCE CRITERIA. 

COMPARE LCB OF LOT WITH LCBFB 

IF THE LCB IS GREATER THAN OR EQUAL TO LCBFB, 

GIVE BONUS TO THE FIRM. 
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According to the algorithm described above, the following are possible events 
for the firms. 

■ Firm may not satisfy our acceptance criteria. This means that firm gets a 
failure vector worse than (1 12 1). 

■ Firm satisfies acceptance criteria, but its LCB may be less than LCBFB. This 
means that firm does not get bonus, but the lot is accepted. 

■ Firm satisfies acceptance criteria, and its LCB may be greater than or equal 
to LCBFB. This means that firm gets the bonus. 

2. Double Sampling Bonus System 

A double sampling plan has an advantage over a single sampling plan. Because 
a double sampling plan involves a larger sample size, it reduces the chance that a man- 
ufacturer who deserves a bonus will not get one. Double Sampling Bonus System per- 
mits the taking of two samples on which to make a decision [Ref 11 : pp. 184-185]. 

In this system, we have two inspection stages. If the firm does not get a bonus 
after the first inspection, then the firm is given a second chance with a second inspection. 
A double sampling bonus system is designated by five numbers. 

«] , «2 » LCBFB, LCB] , LCB2 



where 

■ «, is the sample size for first inspection, ( here «, = 92 ) 

■ is the sample size for second inspection, ( here = 92 ) 

■ LCBFB is the cut off value applied to lower confidence bounds for awarding 

bonuses, 

■ LCBl is the estimated lower confidence bound after first inspection, 

■ LCB2 is the updated estimate of lower confidence bound after second in- 
spection. 
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In the double sampling bonus system, if the firm does not get a bonus after the 
first inspection but does meet the acceptance criteria for the first sample, then a second 
sample is taken. We already have LCB's of possible cases after the first inspection in 
Table 1 1. After the second inspection, we calculate the LCB using an aggregated failure 
vector which includes failures from both samples. When we compute LCB's using ag- 
gregated failure vector after the second inspection, we have a total sample size of 184; 
there are 135 different failure vectors for which the lot meets the acceptance criteria for 
both samples. After tabulating these possible 135 cases and the estimates of R, we used 
the bootstrap procedure to find LCB's for each case. We created 5000 random failure 
vectors for each of them by using case success probabilities in tests. After this, we esti- 
mated MLE's of 5000 failure vectors for each possible case. Finally, we estimated 95 
% LCB of each case. These calculations were done by using programs RAXVEC, 
MLEB, SORT in appendix B, G and C respectively. The results are tabulated at Ap- 
pendix 1. The following algorithm shows us how' double sampling bonus system works 
for pyrotechnic device. 



STEP it 1 
STEP it 2 
STEP it 3 
STEP it 4 
STEP it 5 

STEP it 6 
STEP it 7 

STEP it 8 



DETERMINE LCBFB. 

TAKE A SAMPLE SIZE OF 92 FROM OFFERED LOT. 

APPLY MANUFACTURER AND THREE ENVIRONMENT TESTS. 
COMPARE RESULTS OF TESTS WITH ACCEPTANCE CRITERIA. 
ESTIMATE ITS LOWER CONFIDENCE BOUND ( LCBl ), 

IF THEY MET ACCEPTANCE CRITERIA. 

COMPARE LCBl OF LOT WITH LCBFB 

IF THE LCBl IS GREATER THAN OR EQUAL TO LCBFB, 

GIVE BONUS TO THE FIRM. 

IF THE LCBl IS LESS THAN LCBFB, 
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GIVE A SECOND CHANCE TO THE FIRM FOR BONUS. 



STEP # 9 : 

STEP # 10 : 
STEP if 11 : 

STEP if 12 : 

STEP # 13 : 
STEP if 14 : 



TAKE A NEW SAMPLE SIZE OF 92 FROM OFFERED LOT. 
APPLY MANUFACTURER AND THREE ENVIRONMENT TESTS. 
AFTER GETTING THE NEW FAILURE VECTOR, 

ADD THIS ONE TO THE FIRST FAILURE VECTOR. 
ESTIMATE ITS LOWER CONFIDENCE BOUND ( LCB2 ), 
WITH AGGREGATED FAILURE VECTOR. 

COMPARE LCB2 OF LOT WITH LCBFB 

IF THE LCB2 IS GREATER THAN OR EQUAL TO LCBFB, 

GIVE BONUS TO THE FIRM. 



According to the algorithm described above, the following are possible events 
for the firms. 

■ The firm may not satisfy our acceptance criteria after first inspection. This 
means that firm gets a failure vector worse than (112 1). The lot is not accepted and 
second sample is not taken. 

■ The firm satisfies acceptance criteria, and its LCBl may be greater than or 
equal to LCBFB. This means that it gets bonus after first inspection, and a second 
sample is not needed. 

■ The firm satisfies acceptance criteria, but its LCBl may be less than LCBFB. 
This means that it does not get the bonus after first inspection, but still has a chance to 
get a bonus if it submits a second sample. 

Of those which submit a second sample 

■ LCB2 may be less than LCBFB. It means that firm does not get the bonus, 
but the lot is accepted. 
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■ The aggregated failure vector meets acceptance criteria. But LCB2 may be 
grater than or equal to LCBFB. It means that firm gets bonus. 

3. Multi-Sampling Bonus System 

In this system approach, we may take a sample of size 92 from offered lot three 
times, four times or more. But if we use the three sampling bonus system, after the third 
inspection, we may get(4X4X7x4 = 448 ) 448 different failure vectors. Multi- 
Sampling Bonus System is going to be more computationally intensive with respect to 
double sampling bonus system. This can be done later using with the same reasoning 
in double sampling bonus system. 

C. EXAMPLES 

I. Single Sampling Bonus System 

Let us assume that we decided that lower confidence bound for bonus will be 
0.8500. We are going to use this LCBFB in these examples. 

■ Case # 1 

At the end of inspection, firm has following failure vector 

( 2 0 0 0 ). 

The firm did not meet the acceptance criteria. We immediately reject the lot and no 
bonus is given. 

■ Case # 2 

At the end of inspection, firm has following failure vector 

(10 0 1 ). 

The firm meets the acceptance criteria. We estimate its 95 % lower confidence bound 
to be 
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LCD = 0.8222825 . 



Because LCD < LCBFB , we do not give a bonus to the firm. 

■ Case # 3 

At the end of inspection, firm has following failure vector 

(001 0 j. 

The firm meets the acceptance criteria. We estimate its 95 % lower confidence bound 
to be 



LCB = 0.8766983 . 

Because LCB > LCBFB , we give a bonus to the firm. 

2. Double Sampling Bonus System 
■ Case # 1 

At the end of first inspection, firm has following failure vector. 

(10 0 1 ). 

The firm meets the acceptance criteria. We estimate its 95 % lower confidence bound 
to be 



LCB\ = 0.8222825 . 

Because LCB\ < LCBFB , we give a second chance to the firm. Here is the result of 
second inspection 



( 0 3 2 1 ). 

The aggregated failure vector will be 
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( 1 3 2 2 ). 



It is clear that LCB2 < LCBl for this data, thus we do not give a bonus to the firm. 

■ Case # 2 

At the end of first inspection, firm has following failure vector 

(10 0 1 ). 

The firm meets the acceptance criteria. We estimate its 95 % lower confidence bound 
to be 



LCB\ = 0.8222825 . 

Here, LCB < LCBFB thus we will give a second chance to the firm. Here is the result 
of second inspection 

(001 1 ). 

The aggregated failure vector will be 

(10 1 2 ), 

and we estimate its 95 % lower confidence bound to be 

LCB\ = 0.8369564 . 

Because LCB < LCBFB , we do not give a bonus to the firm. 

■ Case § 3 

At the end of first inspection, firm has following failure vector 

(10 0 1 ). 
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The firm meets the acceptance criteria. We estimate its 95 % lower confidence bound 
to be 



LCB\ = 0.8222825 . 



Because LCB\ < LCBFB , we will give a second chance to the firm. Here is the result 
of second inspection 

( 0 0 0 1 ). 

The aggregated failure vector will be 

( 1 0 0 2 ), 

and we estimate its 95 % lower confidence bound to be 



LCB\ = 0.8706521 . 



Because LCB2 > LCBFB , we will give a bonus to the firm. 
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VI. SIMULATION RESULTS OF BONUS SYSTEMS 



A. BACK GROUND 

We simulated the bonus system to get an idea of how well the bonus system pro- 
posed in the previous chapter works. This provides the user with a means of setting the 
cut off criteria for awarding a bonus. To generate random failure vectors, we need to 
know the probabilities of being successful in each test for the firm. First, we assume that 
the firm has an equal probability of being successful in each test, with the following 
values of P^ = P2 — = P^ 

• 0.9375 

• 0.9500 

• 0.9750 

• 0.9900 

• 0 . 9950 . 

In a second set of simulations, we assume that the manufacturer test has a bigger 
probability of being successful than the other environment tests. The following proba- 
bilities of being successful in manufacturer test were used. 

• 0.9990 

• 0.9950 

• 0.9750 

• 0 . 9500 . 
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In these cases, we used equal probabilities of being successful in other environment tests 
with values between 0.9375 and manufacturer test probability for that given case. 
Finally, we assume that the manufacturer test should have smaller probability of being 
successful. When we study these cases with this assumption, we used probabilities of 
being successful in other environment tests between manufacturer test probability and 
0.9990 for that given case. We assume that minimum value of the manufacturer test 
probability will be 0.9500 because of worst case in last two assumptions. 

We generated 2000 random failure vectors for each possible combination of proba- 
bilities by using program RANVEC in Appendix B. They were used for 1000 repli- 
cations of each bonus system, because the Double Sampling Bonus System ( DSBS ) 
potentially uses two failure vectors per replication. 

After getting the failure vectors, the next step is to decide lower confidence bound 
for giving bonus ( LCBFB ). We chose 0.800, 0.825, 0.850, 0.875, 0.900, 0.950, 0.999 
as LCBFB during our simulations. We simulated bonus systems by using program BO- 
NUS in Appendix G. Program BONUS counts how many times firm gets bonus during 
1000 replications. And then it calculates the bonus percentage dividing counted number 
by 1000. For each scenario this bonus percent is an estimate of the probability of getting 
a bonus. 

B. INITIAL COMPARISON OF SYSTEMS 

In this part, with equal probabilities of being successful in tests, we tried to see the 
difference between the Single and the Double Sampling Bonus Systems. For this reason, 
we simulated Single Sampling Bonus System. Results were tabulated and plotted in next 
pages. 
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Table 12. SINGLE SAMPLING BONUS SYSTEM (EQUAL PROBABILITIES) 



PROB.'S 


CHOSEN LOWER CONFIDENCE BOUNDS FOR BONUS 




0.800 


0.825 


0.850 


0.875 


0.900 


0.950 


0.999 


0.9200 


0.025 


0.007 


0.007 


0.007 


0.002 


0.002 


0.002 


0.9375 


0.060 


0.023 


0.023 


0.023 


0.010 


0.010 


0.007 


0.9450 


0.092 


0.044 


0.044 


0.044 


0.018 


0.018 


0.010 


0.9500 


0.117 


0.053 


0.053 


0.053 


0.022 


0.022 


0.011 


0.9600 


0.222 


0.112 


0.112 


0.112 


0.045 


0.045 


0.028 


0.9700 


0.380 


0.198 


0.198 


0.198 


0.107 


0.107 


0.067 


0.9750 


0.482 


0.281 


0.281 


0.281 


0.163 


0.163 


0.067 


0.9800 


0.609 


0.381 


0.381 


0.381 


0.237 


0.237 


0.106 


0.9850 


0.740 


0.511 


0.511 


0.511 


0.362 


0.362 


0.166 


0.9900 


0.868 


0.646 


0.646 


0.646 


0.502 


0.502 


0.274 


0.9950 


0.955 


0.805 


0.805 


0.805 


0.502 


0.502 


0.407 



It is obvious that there is no difference between some lower confidence bounds for bonus 
from the table above. The firm gets the same bonus percentage when we use 0.825, 
0.850, 0.875 as lower confidence bound for bonus. The same thing occurs when we use 

0. 900 and 0.950 as LCBFB. For this reason, we are going to see four curves in Figure 

1. Figure 2 shows Double Sampling Bonus System with different LCBFB's. 
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SINGLE SAMPUNG BONUS SYSTEM 
DIFFERENT LCB'S FOR BONUS WflTH EQUAL TEST PROBABIUTIES 




Figure 1. Single Sampling Bonus System With Different LCBFB's 



DOUBLE SAMPUNG BONUS SYSTEM 
DIFFERENT LCB'S FOR BONUS WTH EQUAL TEST PROBABIUTIES 




Figure 2. Double Sampling Bonus System With Different LCBFB's 
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As you see in figure 2, we have six difierent curves for Double Sampling Bonus System. 
Because in this system, we can see from the table that LCBFB's 0.900 and 0.950 have 
approximately the same bonus percentage. For this reason we did not plot for 0.950. 



Table 13. DOUBLE SAMPLING BONUS SYSTEM (EQUAL PROBABILITIES) 



PROB.'S 


CHOSEN LOWER CONFIDENCE BOUNDS FOR BONUS 




0.800 


0.825 


0.850 


0.875 


0.900 


0.950 


0.999 


0.9200 


0.026 


0.008 


0.007 


0.007 


0.002 


0.002 


0.002 


0.9375 


0.064 


0.030 


0.025 


0.023 


0.010 


0.010 


0.007 


0.9450 


0.108 


0.059 


0.050 


0.045 


0.020 


0.018 


0.010 


0.9500 


0.142 


0.079 


0.063 


0.057 


0.027 


0.022 


0.011 


0.9600 


0.277 


0.172 


0.144 


0.122 


0.052 


0.045 


0.028 


0.9700 


0.498 


0.351 


0.293 


0.236 


0.138 


0.107 


0.067 


0.9750 


0.624 


0.487 


0.423 


0.353 


0.212 


0.163 


0.106 


0.9800 


0.752 


0.632 


0.565 


0.487 


0.319 


0.237 


0.166 


0.9850 


0.861 


0.799 


0.738 


0.643 


0.486 


0.362 


0.274 


0.9900 


0.948 


0.920 


0.887 


0.799 


0.671 


0.502 


0.407 


0.9950 


0.987 


0.982 


0.970 


0.941 


0.872 


0.695 


0.626 



Figure 2 shows DSBS to be more sensitive in the sense that a higher percentage of 
the firms were awarded a bonus. For this reason, we decided to implement DSBS in all 
simulations. 
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C. SIMULATION RESULTS WITH DIFFERENT LCD'S FOR BONUS 

In this section, results for each chosen lower confidence bound for bonus will be 

presented as follows. We used values in Table 13 to draw plots with equal probabilities 

in each test. Polynomial approximation was used in curve fitting. 

DOUBLE SAMPLING BONUS SYSTEM 
EQUAL TEST PROBABILfTIES: LCBEB -= O.BOO 




Figure 3. Double Sampling Bonus System With LCBFB = 0.800 



Table 14. DSBS (EQUAL PROBABILITIES) LCBFB = 0.800 



LOWER CONFIDENCE BOUND FOR BONUS IS 0.800 


PROB.'S 


BONUS % 


PROB.'S 


BONUS % 


PROB.'S 


BONUS % 


0.9200 


0.026 


0.9375 


0.064 


0.9450 


0.108 


0.9500 


0.142 


0.9600 


0.277 


0.9700 


0.498 


0.9750 


0.624 


0.9800 


0.752 


0.9850 


0.861 


0.9900 


0.948 


0.9950 


0.987 
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In Figure 3, when the probabilities of being successful in each test increase, then the 
bonus percent increases. For example, when the test probabilities is equal to 0.9200, 
then the bonus percent is 0.026. If the firm increases its probabilities of being succesful 
in each test to 0.9900, then the bonus percent becomes 0.948. The result of the double 
sampling bonus system with different lower confidence bounds are tabulated and plotted 
in Appendix J. 

To see whether the probability of passing the manufacturer test effects the bonus 
percentage differently than the other probabilities, we assumed that the environment 
tests would have equal probabilities. But the firm will have a different probability of 
being successful in the manufacturer test. In the following table, the first row represents 
probability of being successful in the environment tests and the first column represents 
probability of being successful in the manufacturer test. The intersection of rows and 
columns gives us the bonus percentage of a firm with given probabilities of being suc- 
cessful in the tests. These procedures were done for each LCBFB value separately. The 
results are tabulated and plotted in Appendix K. Bonus percentages for LCBFB = 
0.8000 are plotted on the following page. It is clear that the bonus percentage of the 
firm will be high if the firm has big probabilities of being succesful in both manufacturer 
test and other joint environment tests. 
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Bonus percentages are tabulated and ploted with difierent probabilities. 



LCBFB = 0.800 



DOUBLE SAMPLING BONUS SYSTEM 
LCBFB = 0.800 




Figure 4. Double Sampling Bonus System With LCBFB = 0.800 



Table 15. DSBS (DIFFERENT PROBABILITIES) LCBFB = 0.800 





0.950 


0.975 


0.990 


0.995 


0.950 


0.142 


0.476 


0.683 


0.700 


0.975 


0.176 


0.624 


0.869 


0.891 


0.990 


0.202 


0.688 


0.948 


0.973 


0.995 


0.206 


0.698 


0.962 


0.987 
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D. BONUS PERCENTAGE ( BPRCT ) FORMULATION 



In this section we tty to approximate probability of getting a bonus as a function 
of lower confidence bound for bonus ( LCBFB ), probability of passing manufacturer 
test ( P| ) , and probability of passing environment tests ( P 2 )• To do this, we use re- 
gression analysis with GRAFSTAT using the bonus percentages which are the simulated 
values of the probability of getting a bonus. If a reasonable relationship between the 
bonus percentage and LCBFB and the probabilities P, and P-^ can be found; it can be 
used to set LCBFB without resorting to simulation. 

We want to formulate bonus percentages as a function of the following: 

■ Lower confidence bound for bonus ( LCBFB ) 

■ Probability of passing manufacturer test P, 

■ Probability of passing environment tests Pj. 

After polynomial approximation, we can see that plotted graphs (Figures 3-10) look 
like logistic growth curves (Ref 12: p. 383]. Then we can formulate bonus probabilities 
(BP) using the logistic growth function. 



BP 



1 

1 I ^ 

1 -b e 



(7.4) 



where 



/I = /?0 + /?, P, + /?2 P '1 + fi^P2 + Pa pI + Ps Px P2 + P6 LCBFB (7.5) 
We can make a transformation described as below. 



1 -f e 



A 



1 

BP 



(7.6) 



A _ 1 - BP 

^ BP 



(7.7) 
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(7.8) 



Now we have a linear equation as a function of LCBFB, f, and P^. We may use 
simulation results and approximate parameters doing a linear regression. 

■ Let Y be the ( n x 1) column vector of observations on dependent variable. 

■ Let X be the ( n x /j' ) matrix consisting of a column of ones, which is labeled 1, 
followed by p column vectors of observations on independent variables. 

■ Let P be the ( /?' x 1 ) vector of parameters to be estimated. 

■ Let c be the ( n x 1 ) vector of random errors. Then 

Y = Xp + c (7.9) 

We can obtain observations on dependent variable as below, 

V , r 1 - BPRCT 1 

^ L BPRCT J 

where BPRCT is the bonus percentage obtained from the simulation. As an example 
we have the following information from simulations. 

. LCBFB = 0.875 
. P, = 0.9500 
. P\ = 0.9025 
. Pj = 0.9375 
. PI = 0.8789 
. P, Pj = 0.8910 
. BPRCT = 0.029 
. Y = 3.5110. 

We can write following equation. 
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3.5110 = Po + jff, 0.95 + /?2 0.9025 4- ^3 0.9375 + /?4 0.8789 + /?5 0.8910 4- /?6 0.875 + £. 



Wc deliberately chose 30 random results from simulations. We formulate them as 
above. We can estimate unknown parameters on these equations doing linear re- 
gression. 30 equations can be written with matri.x notation as follows: 



y = 



rO) 

V(2) 

V(3) 



V(30) 



(7.11) 



And X matrix is designed as follows: 



^.( 1 ) 


^?( 1 ) 




Plii) 


P, P2(i) 




^i( 2 ) 


Px{2) 


P2(2) 


P^(2) 


Pi P2(2) 




P,(3) 


Pxi^) 


P 2 O) 


P^(V 


Pi P 2 O) 


(7.12) 


/’,(30) 


/>?(30) 


P^{30) 


P^(30) 


Pi ^2(30) 





And parameter vector will be as below. 



P = 



Po 

Px 

P2 

P3 

P, 

Ps 

P6 



(7.13) 
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We estimate parameters with linear regression using GRAFSTAT packages. And 
then we formulate the bonus percentage using the estimates of parameters. The esti- 
mated bonus probabilities are given by 



BP 



1 



1 + e 






(7.14) 



The following estimates of the parameters were obtained from the regression. 



Table 16. COEFFICIENT VECTORS OF REGRESSION ANALYSIS 



PARAMETERS 


ESTIMATED VALUE 


Po 


-1042.3 


p. 


695.11 


p, 


-86.661 


p. 


1560 


ps 


-569.81 


p. 


-570.62 


Pi 


12.266 



The standard error was 0.45471 after the linear regression. We can use these esti- 
mates of parameters in equation (7.5). And we can estimate the bonus probability of the 
firm when the probability of passing tests and lower confidence for bonus are known. 
Some examples are calculated with equation (7.5). Results of the calculations and 
comparison with simulation are summarized in following table. 
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Table 17. BONUS PERCENTAGES WITH REGRESSION ANALYSIS 



LCBFB 


A 


Pz 


SIMU- 

LATION 


EST.VALUE 


0.800 


0.9750 


0.9750 


0.624 


0.564 


0.800 


0.9990 


0.9375 


0.106 


0.062 


0.825 


0.9750 


0.9950 


0.882 


0.911 


0.825 


0.8500 


0.9750 


0.064 


0.077 


0.850 


0.9500 


0.9900 


0.612 


0.615 


0.850 


0.9900 


0.9750 


0.476 


0.529 


0.875 


0.9000 


0.9750 


0.114 


0.080 


0.875 


0.9700 


0.9700 


0.236 


0.211 


0.900 


0.9950 


0.9990 


0.985 


0.945 


0.900 


0.8750 


0.9900 


0.158 


0.096 


0.950 


0.9450 


0.9450 


0.018 


0.010 


0.950 


0.9750 


0.9990 


0.839 


0.790 


0.999 


0.9850 


0.9850 


0.274 


0.335 


0.999 


0.9900 


0.9250 


0.004 


0.002 
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VII. CONCLUSIONS AND RECOMMENDATIONS 



We estimated the reliability of pyrotechnic device from the sampling plan using the 
following models: 

■ Maximum likelihood assuming independence of pyrotechnic device activation in 
different environments; 

■ Log linear model incorparating some dependence between pyrotechnic device ac- 
tivation in different environments; 

■ Worst case scenario which gives a lower bound for the estimated reliability for the 
general model where no assumptions are made about the form of dependence between 
pyrotechnic device activation in different environments. 

Using these models and based on sampling plan data, estimates of overall reliability 
along with 95 % lower confidence bound are obtained. We computed the lower confi- 
dence bound for each possible case by bootstrapping. 

Results from the first model are not consistent with the way pyrotechnic devices 
operate. In particular, the estimated reliabilities are not ordered as we expected them 
to be. Intiutively, we expect samples with fewer failures to give smaller R values and 
corresponding LCD's than samples with more failures. This is not the case for the first 
model. For example, the failure vector that has the maximum number of failures for 

A 

each test ( a total of 5 failures ) is in the middle of the order with respect to R and the 
95 % lower confidence bound. The failure vector (0101) with a total of 2 failures 

A 

has a much lower R and LCD. The discrepency between the results of this model and 
what we expected to see are probably due to the fact that there is a dependence between 
the events which a device activates under different environments. 
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The results of the log linear model design with two-way dependence assumptions are 
similar to the first model. Two-way interaction terms were used in this model. More- 
over, reliabilities of several cases were not calculated because of computational limita- 
tions. Due to the extensive censoring in the sample data, it does not appear to be 
possible to estimate the reliability based on models which incorparate dependence. 
Thus, we turn to finding lower bounds for the estimated reliability based on models with 
dependence. 

Finally, the worst case scenario model gives the most reasonable results for this 
problem. Both the lower bounds for estimates of reliabilities and 95 % lower confidence 
bounds for these lower bounds are ordered according to the total number of failures. 
Thus, the results from the worst case scenario were used to implement the bonus system 
for pyrotechnic devices. 

After getting the estimate of reliability and 95 % lower confidence bound for each 
case, we tried to design a bonus system to improve the quality of pyrotechnic devices. 
We used 95 % lower confidence bounds instead of the estimated overall reliabilities to 
decide whether to give bonuses. Two sampling plans for giving the bonus to the man- 
ufacturers were considered: 

■ Single Sampling Bonus System 

■ Double Sampling Bonus System. 

We simulated two sampling bonus systems to see the difference between them. 
We concluded that the double sampling system is more sensitive than single sampling 
bonus system. 

To formulate an approximate bonus percentage as a function of lower confidence 
bound for bonus ( LCBFB ), the probability of passing manufacturer test P^, and the 
probability of passing the environment tests P^, we used regression analysis with 
GRAFSTAT. 
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Because the sampling plan results in so much censoring, the only reasonable esti- 
mates of overall reliability that we obtained were actually lower bounds. This makes the 
bonus system conservative in the sense that a bonus might not be awarded what it is 
deserved. Therefore, if a bonus system is to be implemented, a more comprehensive 
sampling plan needs to be devised which allows estimation of R. A simple solution to 
this problem can be to apply all environmental tests to the same sample which would 
give a measure of dependence between environmental tests. 



76 



APPENDIX A. PROGRAM MLEA 




THIS IS A FORTRAN PROGRAM TO CALCULATE THE RELIABILITY OF AN 
ITEM AFTER EXPOSURE TO SEVERAL ENVIRONMENTS WITH INDEPENDENCE 
ASSUMPTION WHICH IS DESCRIBED IN CHAPTER I . THE PROGRAM READS 
5000 SUCCESS VECTORS, WHICH ARE RANDOMLY GENERATED BY PROGRAM 
RANVEC, FROM AN INPUT FILE CALLED SUCVECT ONE BY ONE. AFTER 
CALCULATION IT THEN PROGRAM WRITES ESTIMATED RELIABILITIES TO 
AN OUTPUT FILE CALLED RESULT. 




SOM : NUMBER OF SUCCESSFUL ITEMS IN MANUFACTURER TEST. 

SOTH : NUMBER OF SUCCESFUL ITEMS IN TEMPERATURE AND HUMIDITY 
TEST. 

SOV ; NUMBER OF SUCCESSFUL ITEMS IN VIBRATION TEST. 

SOA ; NUMBER OF SUCCESSFUL ITEMS IN ALTITUDE TEST. 

RIH : ESTIMATED PROBABILITY OF PASSING FROM MANUFACTURER TEST. 

R2H : ESTIMATED PROBABILITY OF PASSING FROM TEMPERATURE AND 

HUMIDITY TEST. 

R3H : ESTIMATED PROBABILITY OF PASSING FROM VIBRATION TEST 

R4H : ESTIMATED PROBABILITY OF PASSING FROM ALTITUDE TEST 

RHMLE : ESTIMATED RELIABILITY OF ITEM AFTER EXPOSURE TO 
SEVERAL ENVIRONMENT TESTS. 

N : SAMPLE SIZES FOR EACH TEST 
X : NUMBER OF SUCCESSFUL ITEMS IN EACH TEST. 

FLAG : INDICATOR VARIABLE FOR DETERMINING EASY AND HARD CASE. 
RIMAX : RIH VALUE WHICH MAXIMIZES LIKELIHOOD. 



* 



* 



* 



R2MAX : R2H VALUE WHICH MAXIMIZES LIKELIHOOD. 

R3MAX : R3H VALUE WHICH MAXIMIZES LIKELIHOOD. 

R4MAX : R4H VALUE WHICH MAXIMIZES LIKELIHOOD. 

LMAX : HIGHEST MAXIMUM LIKELIHOOD VALUE. 

L : LIKELIHOOD VALUES AT THE END OF EACH HARD CASE. 

V : TOTAL SAMPLE SIZES IN EACH HARD CASE. 

QIH : REORGANIZED PROBABILITY OF MANUFACTURER TEST. 

Q2H : REORGANIZED PROBABILITY OF TEMP. AND HUMIDITY TEST. 

Q3H : REORGANIZED PROBABILITY OF VIBRATION TEST. 

Q4H : REORGANIZED PROBABILITY OF ALTITUDE TEST. 



* 



TYPE DECLARATION 

REAL SOM(5000),SOTH(5000),SOV(5000),SOA(5000),R1H,R2H,R3H,R4H, 
+RHMLE(5000),N(4),X(4),R1MAX,R2MAX,R3MAX,R4MAX, LMAX, V(7), QIH, Q2H, 
+Q3H,Q4H,NTOT,L(7) 

INTEGER I,J 
LOGICAL FLAG(4) 






READING SUCCESS VECTORS FROM SUSVECT FILE 
DO 60 1=1,5000 

READ(7,*) SOM(I),SOTH (I) ,SOV(I) ,SOA(I) 

■ A>W. - A -- /. - vV - Av\ - .V/. ' A^VA - J. ' .VA ' <VAAAA Tt * - ** - A ' ^**A Tfr A/. - y.»V Vt A ' A - vVvV%^^ ^ ^ ^ 

FILES FOR READING AND WRITING 

CALL EXCMS (’FILEDEF 7 DISK SUCVECT DATA Al') 

CALL EXCMS ('FILEDEF 16 DISK RESULT DATA Al’) 

A A A A A A A T f. ' ^. - A A -Ar A-A - A - A - A T. V A' A 'A- A A -Ari' f A - A^.^V A ' A - A-A 

INITIALIZATION OF SAMPLE SIZES 

N(l)= 20.0 

N(2)= 20.0 

N(3)= 32.0 

N(4)= 20.0 

INITIALIZATION OF FLAG VARIABLES 
DO 10 J=2,4 

FLAG(J)= .FALSE. 
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10 



* 



* 

* 



CONTINUE 

NUMBER OF SUCCESSES IN EACH TEST 
X(l)= SOM(I) 

X(2)= SOTH(I) 

X(3)= SOV(I) 

X(4)= SOA(I) 



CHECK OPERATION FOR EASY CASE 
IF ((X(l)/20.00).GE.(X(2)/20.00)) THEN 
FLAG(2)= .TRUE. 

END IF 

IF ((X(l)/20.00).GE.(X(3)/32.00)) THEN 
FLAG(3)= .TRUE. 

END IF 

IF ((X(l)/20.00).GE.(X(4)/20.00)) THEN 
FLAG(4)= .TRUE. 

END IF 

IF ( FLAG(2).AND.FLAG(3).AND.FLAG(4) ) THEN 



ir.' c is* *irlrie*it ***** - J. - ****M. - * -k*ii 



*********************************** 



CALCULATIONS IN EASY CASE 
R1H= X(l)/20.0 
R2H= X(2)/20.0 
R3H= X(3)/32.0 
R4H= X(4)/20.0 
Q1H= RIH 
Q2H= R2H/Q1H 
Q3H= R3H/Q1H 
Q4H= R4H/Q1H 

RHMLE(I)= Q1H*Q2H*Q3H*Q4H 
GO TO 50 
END IF 



CALCULATIONS IN HARD CASES 
INITIALIZATION 
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DO 30 J=l,7 
V(J)= 0.0 
L(J)= 0.0 
30 CONTINUE 
LMAX= 0.0 
NTOT= 0.0 



40 



* 



CASE 1 RIH.LT.1.0 AND RIH = R2H = R3H = RAH 
DO 40 J=l,4 

V(l)= V(1)+X(J) 

NTOT= NTOT+N(J) 

CONTINUE 
R1H= V(l)/NTOT 

L( 1)= ( (R1H)**V( 1) )*( (NTOT-V( 1) )/NTOT)**(NTOT-V( 1) ) 
IF (L(l).GT.LMAX) THEN 
LMAX= L(l) 

R1MAX= RIH 
R2MAX= RIH 
R3MAX= RIH 
R4MAX= RIH 
END IF 



*** 



CASE 2 R1H.lt. 1.0, RIH = R2H = R3H , RAH IS BETWEEN 0.0 AND RIH 
IF ( X(4)/20.0.NE.1.0 ) THEN 
V(2)= X(l)+X(2)+X(3) 

NTOT= N(l)+N(2)+N(3) 

R1H= V(2)/NTOT 

IF((X(4)/20.00).LE.R1H) THEN 

L(2)= ((RlH)**V(2))*((NTOT-V(2))/NTOT)**(NTOT-V(2))*((X(4) 
+/20.0)**X(4))*(((20.0-X(4))/20.0)**(20.0-X(4))) 

IF (L(2).GT.LMAX) THEN 



LMAX= L(2) 
R1MAX= RIH 
R2MAX= RIH 
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R3MAX= RIH 
R4MAX= X(4)/20.0 
END IF 
END IF 
END IF 



* 

* 



CASE 3 RIH. LT. 1.0, RIH = R2H = R4H , R3H IS BETWEEN 0.0 AND RIH 
IF ( X(3)/32.0.NE.1.0 ) THEN 
V(3)= X(l)+X(2)+X(4) 

NTOT= N(l)+N(2)+N(4) 

R1H= V(3)/NTOT 

IF ((X(3)/32.00).LE.R1H) THEN 

L(3)= ((RlH)**V(3))*((NTOT-V(3))/NTOT)**(NTOT-V(3))*((X(3) 
+/32.0)**X(3))*(((32.0-X(3))/32.0)**(32.0-X(3))) 

IF (L(3).GT.LMAX) THEN 
LMAX= L(3) 

R1MAX= RIH 
R2MAX= RIH 
R3MAX= X(3)/32.0 
R4MAX= RIH 
END IF 
END IF 
END IF 

CASE 4 RIH. LT. 1.0, RIH = R3H = R4H , R2H IS BETWEEN 0.0 AND RIH 
IF ( X(2)/20.0.NE.1.0 ) THEN 
V(4)= X(l)+X(3)+X(4) 

NTOT= N(l)+N(3)+N(4) 

R1H= V(4)/NT0T 

IF ((X(2)/20.00).LE.R1H) THEN 

L(4)= ( (R1H)**V(4) )*( (NT0T-V(4) )/NT0T)**(NT0T-V(4) )*( (X(2) 
+/20.0)**x(2))*(((20.0-X(2))/20.0)**(20.0-X(2))) 

IF (L(4).GT.LMAX) THEN 
LMAX= L(4) 
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R1MAX= RIH 
R2MAX= X(2)/20.0 
R3MAX= RIH 
R4MAX= RIH 
END IF 
END IF 
END IF 

CASE 5 RIH. LT. 1.0, RIH = R2H R3H , R4H ARE BETWEEN 0.0 AND RIH 
IF ((X(3)/32.0.NE.1.0).AND.(X(4)/20.0.NE.1.0)) THEN 
V(5)= X(l)+X(2) 

NTOT= N(l)+N(2) 

R1H= V(5)/NTOT 

IF((X(3)/32.00) .LE.R1H.AND.(X(4)/20.00).LE.R1H) THEN 

L(5)= ((R1H)**V(5))*((NT0T-V(5))/NT0T)**(NT0T-V(5))*((X(3) 
+/32.0)**X(3))*((32.0-X(3))/32.0)**(32.0-X(3))*((X(4)/20.0)** 
+x(4))*((20.0-X(4))/20.0)**(20.0-X(4)) 

IF (L(5).GT.LMAX) THEN 
LMAX= L(5) 

R1MAX= RIH 
R2MAX= RIH 
R3MAX= X(3)/32.0 
R4MAX= X(4)/20.0 
END IF 
END IF 
END IF 

CASE 6 RIH. LT. 1 . 0 , RIH = R3H R2H , R4H ARE BETWEEN 0.0 AND RIH 
IF ((X(2)/20.0.NE.1.0).AND.(X(4)/20.0.NE.1.0)) THEN 
V(6)= X(l)+X(3) 

NTOT= N(l)+N(3) 

R1H= V(6)/NTOT 

IF((X(2)/20.00) .LE.R1H.AND.(X(4)/20.00).LE.R1H) THEN 
L(6)= ((R1H)**V(6))*((NT0T-V(6))/NT0T)**(NT0T-V(6))*((X(2) 



+/20.0)**X(2))*((20.0-X(2))/20.0)**(20.0-X(2))*((X(4)/20.0)- 
+x(4))*((20.0-X(4))/20.0)**(20.0-X(4)) 

IF (L(6).GT.LMAX) THEN 
LMAX= L(6) 

R1MAX= RIH 
R2MAX= X(2)/20.0 
R3MAX= RIH 
R4MAX= X(4)/20.0 
END IF 
END IF 
END IF 

CASE 7 RIH. LT. 1.0, RIH = R4H R2H , R3H ARE BETWEEN 0.0 AND RIH 
IF ((X(2)/20.0.NE.1.0).AND.(X(3)/32.0.NE.1.0)) THEN 
V(7)= X(l)+X(4) 

NTOT= N(l)+N(4) 

R1H= V(7)/NTOT 

IF((X(2)/20.00).LE.RlH.AND.(X(3)/32.00) .LE.RIH) THEN 
L(7)= ((RlH)**V(7))*((NTOT-V(7))/NTOT)**(NTOT-V(7))*((X(3) 
+/32.0)**X(3))*((32.0-X(3))/32.0)**(32.0-X(3))*((X(2)/20.0)**X(2)) 
+*((20.0-X(2))/20.0)**(20.0-X(2)) 



IF (L(7) .GT.LMAX) THEN 
LMAX= L(7) 

R1MAX= RIH 
R2MAX= X(2)/20.0 
R3MAX= X(3)/32.0 
R4MAX= RIH 
END IF 
END IF 
END IF 
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CALCULATION OF RHMLE VALUE WITH RESPECT TO CASE WHICH HAS LARGEST 

MAXIMUM LIKELIHOOD IN HARD CASE 

R1H= RIMAX 

R2H= R2MAX 

R3H= R3MAX 

R4H= R4MAX 

Q1H= RIH 

Q2H= R2H/Q1H 

Q3H= R3H/Q1H 

Q4H= R4H/Q1H 

RHMLE(I)= Q1H*Q2H*Q3H*Q4H 















J I — 2 V J1 1 t U. 



WRITING AFTER EACH CALCULATION 
50 WRITE (16,*) RHMLE(I) 

60 CONTINUE 
STOP 
END 
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APPENDIX B. PROGRAM RANVEC 



PROGRAM RANVEC 

THIS IS THE PROGRAM TO GENARETE RANDOM NUMBERS FROM BINOMIAL 
DISTRIBUTION. THE PROGRAM READS PROBABILITIES OF BEING SUCCESSFUL 
IN FOUR TESTS INTERACTIVELY. IT GENERATES UNIFORMLY DISTRIBUTED 
RANDOM NUMBERS WITH THESE PROBABILITIES ACORDING TO SAMPLE SIZE 
OF EACH TEST. FOR EACH TEST PROGRAM COUNTS UNIFORMLY DISTRIBUTED 
RANDOM NUMBERS, WHICH HAVE GREATER THAN OR EQUAL PROBABILITY 
WITH RESPECT TO THE GIVEN PROBABILITY FOR THAT TEST. TOTAL COUNTS 
GIVE US SUCCESSFUL ITEM NUMBERS FOR EACH TEST. PROGRAM UPDATES 
SEEDS AND CALLS SUBROUTINE RANNUM IN EACH ITERATION. THE PROGRAM 
GENERATES 5000 SUCCESS VECTOR AND WRITES THEM TO AN OUTPUT FILE 
CALLED SUSVECT. 

VARIABLES • 



: PROBABILITY OF SUCCESS IN MANUFACTURER TEST. 

: PROBABILITY OF SUCCESS IN TEMP. AND HUMIDITY TEST. 

: PROBABILITY OF SUCCESS IN VIBRATION TEST 

: PROBABILITY OF SUCCESS IN ALTITUDE TEST 

: COUNTER FOR MANUFACTURER TEST 
: COUNTER FOR TEMP. AND HUMIDITY TEST 
: COUNTER FOR VIBRATION TEST 
: COUNTER FOR ALTIUTDE TEST 

: BINOMIALLY DISTRIBUTED RANDOM NUMBER FOR MANUFACTURER 
TEST 

: BINOMIALLY DISTRIBUTED RANDOM NUMBER FOR TEMP. AND 
HUMUDITY TEST 



PSIM 

PSITH 

PSIV 

PSIA 

NUMl 

NUM2 

NUM3 

NUM4 

X 

V 



* 



* 



* 



y : BINOMIALLY DISTRIBUTED RANDOM NUMBER FOR VIBRATION TEST 

Z : BINOMIALLY DISTRIBUTED RANDOM NUMBER FOR ALTITUDE TEST 

A : UNIFORMLY DISTRIBUTED RANDOM NUMBER FOR MANUFACTURER 

TEST 

B : UNIFORMLY DISTRIBUTED RANDOM NUMBER FOR TEMP. AND 

HUMIDITY TEST 

C : UNIFORMLY DISTRIBUTED RANDOM NUMBER FOR VIBRATION TEST 

D : UNIFORMLY DISTRIBUTED RANDOM NUMBER FOR ALTITUDE TEST 

ISEED : SEED NUMBER FOR MANUFACTURER TEST 

KSEED : SEED NUMBER FOR TEMP. AND HUMIDITY TEST 

LSEED : SEED NUMBER FOR VIBRATION TEST 

MSEED : SEED NUMBER FOR ALTITUDE TEST 



TYPE DECLARATION 

REAL PSIM,PSITH,PSIV,PSIA,NUM1,NUM2,NUM3,NUM4,X(5000) ,V(5000) 
+,Y(5000) ,Z(5000) ,A,B,C,D 
INTEGER ISEED, KSEED, LSEED, MSEED 



* INITIALIZATION 
ISEED = 45267 
KSEED = 113234 
LSEED = 435 
MSEED = 1 

* READING TEST PROBABILITIES OF BEING SUCCESSFUL IN EACH TEST 
WRITE(*,*) 'PLEASE WRITE PROBABILITY OF BEING SUCCESSFUL IN 

+ MANUFACTURER TEST* 

READ (*,*) PSIM 

WRITE(*,*) 'PLEASE WRITE PROBABILITY OF BEING SUCCESSFUL IN 
+ TEMPERATURE AND HUMIDITY TEST’ 

READ (*,*) PSITH 

WRITE(*,*) 'PLEASE WRITE PROBABILITY OF BEING SUCCESSFUL IN 
+ VIBRATION TEST' 

READ (*,*) PSIV 
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WRITE(*,*) 'PLEASE WRITE PROBABILITY OF BEING SUCCESSFUL IN 
+ ALTITIDE TEST' 

READ (*,*) PSIA 



OPENING AN OUTPUT FILE TO WRITE THE RESULTS 
CALL EXCMS ( 'FILEDEF 16 DISK SUSVECT DATA Al') 



GENERATION 
DO 50 J=l, 5000,1 






INITIALIZATION IN EACH ITERATION 
NUM1= 0.0 
NUM2= 0.0 
NUM3= 0.0 
NUM4= 0.0 



SEEDS UPDATATION IN EACH ITERATION 
ISEED=ISEED+17 
KSEED=KSEED+1356 
LSEED=LSEED+1 
MSEED=MSEED+789 

GENERATION OF SUCCESSFUL NUMBER OF ITEMS IN MANUFACTURER TEST 
DO 10 1=1,20,1 

CALL RANNUM ( 1 , ISEED,0 . 0 , 1 . 0 ,0 . 0 , A) 

IF ( A.LT.PSIM ) THEN 
NUM1=NUM1+1 
END IF 
X(J) = NUMl 
10 CONTINUE 



GENERATION OF SUCCESSFUL NUMBER OF ITEMS IN TEMP. AND HUM. TEST 
DO 20 K=l,20,l 

CALL RANNUM ( 1 ,KSEED,0 . 0 , 1 . 0, 0 . 0 ,B) 



87 



IF ( B.LT.PSITH ) THEN 
NUM2=NUM2+1 
END IF 
V(J) = NUM2 
20 CONTINUE 



GENERATION OF SUCCESSFUL NUMBER OF ITEMS IN VIBRATION TEST 
DO 30 L=l,32,l 

CALL RANNUM ( 1 ,LSEED, 0 . 0 , 1 . 0 , 0 . 0 ,C) 

IF ( C.LT.PSIV ) THEN 
NUM3=NUM3+1 
END IF 
Y(J) = NUM3 
30 CONTINUE 



GENERATION OF SUCCESSFUL NUMBER OF ITEMS IN ALTITUDE TEST 
DO 40 M=l,20,l 

CALL RANNUM ( 1 ,MSEED, 0 . 0 , 1 . 0 ,0 . 0 ,D) 

IF ( D.LT.PSIA ) THEN 
NUM4=NUM4+1 
END IF 
Z(J) = NUM4 
40 CONTINUE 



WRITING THE RESULTS TO AN OUTPUT FILE AS 4 TUPLE 
WRITE (16,1) X(J),V(J),Y(J),Z(J) 

1 FORMAT (1X,F12.7,4X,F12.7,4X,F12.7,4X,F12.7,4X) 

50 CONTINUE 
STOP 
END 



* 



SUBROUTINE RANNUM(DISTN, SEED, RPARMl, RPARM2, IPARM, X) 


















THIS SUBROUTINE IS A PART OF SIMUTIL FORTRAN WHICH IS WRITTEN 
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* 



* 



* 



* 



•k 



k 



k 



BY DR. M. P. BAILEY. THIS SUBROUTINE PROVIDES AN INTERFACE WITH 
THE LLRANDOMII ROUTINES PROVIDED IN THE NONIMSL LIBRARY. THE 
PARAMETER REQUIRMENTS AND CALLING PROCEDURES ARE AS FOLLOWS: 
DISTN = DISTRIBUTION TYPE YOU WANT TO SELECT AN INTEGER BETWEEN 1 
AND 7. 

SEED = THE RANDOM NUMBER SEED YOU WISH TO USE. 

RPARMl, RPARM2, AND IPARM ARE REAL AND INTEGER PARAMETERS PASSED 
TO THE ROUTINE WITH MEANINGS WHICH VARY WITH THE TYPE OF DISTRI_ 
BUTION YOU DESIRE. 

X = THE RETURNED RANDOM NUMBER, IT IS ALWAYS REAL. 

DISTRIBUTION NUMBERS AND THE ASSOCIATED PARM DEFINITIONS 

1- -UNIFORM ON THE INTERVAL RPARMl TO RPARM2. 

2- -NORMAL WITH MEAN RPARMl AND VARIANCE RPARM2. 

3- -EXPONENTIAL WITH RATE RPARMl. 

4- -C0UCHY WITH A = RPARMl AND B = RPARM2. 

5- -GAMMA WITH SHAPE RPARM2 AND RATE RPARMl. 

6- -POISSON WITH RATE RPARMl. 

7- -GEOMETRIC WITH P = RPARMl. 



TYPE DECLARATION 

REAL RPARMl ,RPARM2 ,X , TEMP , VARIAT( 1 ) 

INTEGER DISTN, SEED, IPARM, N 

k k kir kkk kk k kkkkkkkkkk kr k k- .Ukifkk kk kk k kkkkkkkkkkk kkkk k -k kir kkk k kk kkk -kkk 

IF (DISTN. LE.O. OR. DISTN. GT. 8) THEN 

WRITEdO, *) 'ILLEGAL CALL TO RANNUM, BAD DISTN' 

STOP 

ENDIF 

GOTO (10, 20, 30, 40, 50, 60, 70), DISTN 



GENERATE A UNIFORM BETWEEN RPARMl AND RPARM2 
10 CONTINUE 

IF (RPARMl - RPARM2.EQ.0) THEN 

WRITEdO, *) 'ILLEGAL EQUAL RPARMS IN RANNUM' 
STOP 
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20 



ENDIF 

IF (RPARM1.GT.RPARM2) THEN 
TEMP = RPARMl 
RPARMl = RPARM2 
RPARM2 = TEMP 
ENDIF 

CALL LRND(SEED, VARIAT, 1, 1, 0) 

VARIAT(l) = RPARMl + (RPARM2 - RPARMl) * VARIAT(l) 



GOTO 80 









GENERATE A NORMAL WITH MEAN RPARMl AND STDDEV RPARM2 
CALL LNORM(SEED, VARIAT, 1, 1, 0) 

VARIAT(l) = (VARIAT(l) * RPARM2) + RPARMl 



GOTO 80 












* GENERATE AN EXPONENTIAL WITH RATE (1/MEAN) RPARMl 
30 CONTINUE 

IF (RPARMl. EQ.O) THEN 

WRITEdO, *) 'ILLEGAL ZERO RATE IN RANNUM' 

STOP 

ENDIF 

CALL LEXPN(SEED, VARIAT, 1, 1, 0) 

VARIAT(l) = VARIAT(l) / RPARMl 
GOTO 80 

* GENERATE A COUCHY WITH A = RPARMl AND B = RPARM2 
40 CONTINUE 

IF (RPARM2.LE.0) THEN 

WRITEdO, *) 'ILLEGAL COUCHY SPREAD IN RANNUM, B = ' ,RPARM2 
STOP 
ENDIF 

CALL LCCHY(SEED, VARIAT, 1, 1, 0) 

VARIAT(l) = (VARIAT(l) * RPARM2) + RPARMl 
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GOTO 80 



* 



50 



* 



GENERATE GAMMA WITH SHAPE RPARM2 AND RATE RPRAMl 
CONTINUE 

IF (RPARMl.LE.O) THEN 

WRITEdO, *) 'ILLEGAL NONPOSITIVE GAMMA RATE IN RANNUM' 
STOP 
ENDIF 

IF (RPARM2.LE.0) THEN 

WRITEdO, *) 'ILLEGAL SHAPE PARAMETER IN RANNUM' 

STOP 

ENDIF 

CALL LGAMA(SEED, VARIAT, 1, 1, 0, RPARM2) 

VARIAT(l) = VARIAT(l) * (1.0 / RPARMl) 

GOTO 80 









60 



GENERATE POISSON WITH RATE RPRAMl 
CONTINUE 

IF (RPARMl.LE.O) THEN 

WRITEdO, *) 'ILLEGAL POISSON RATE IN RANNUM' 
STOP 
ENDIF 

CALL LPOIS(SEED, VARIAT, 1, 1, 0, RPARMl) 

GOTO 80 



GENERATE GEOMETRIC WITH P = RPRAMl 
70 CONTINUE 

IF (RPARMl.LE.O) THEN 

WRITEdO, *) 'ILLEGAL GEOM PROB IN RANNUM' 
STOP 
ENDIF 

CALL LGEOM(SEED, VARIAT, 1, 1, 0, RPARMl) 
GOTO 80 
80 CONTINUE 
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X = VARIAT(l) 
END 
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APPENDIX C. PROGRAM SORT 



* 



* 



* 



PROGRAM SORT 

-il^ ^ r b_^ t ^ ^ ^ ^ * * ^ * ^ ^ ^ ^ mS^ mS^ mX^ mi ^ mSm mJ^^im mJ^^immX^ 

THIS IS THE SORTING PROGRAM. PROGRAM USES BUBBLE SORT ALGORITHM. 
PROGRAM READS ESTIMATED RELIABILITIES FROM AN INPUT FILE CALLED 
RESULT. IT SORTS FROM SMALLEST TO LARGEST, AND WRITES IN TO AN 
OUTPUT FILE CALLED FRESULT WITH 95 % LOWER CONFIDENCE BOUND. 

VARIABLES 



A : ESTIMATED RELIABILITY 

FLAG : INDICATOR VARIABLE TAKES VALUE ' OF ' AND ’ OFF ’ 



TYPE DECLARATION 
CHARACTER FLAG*3 
REAL A(5000) 
INTEGER I,N,J 



OPENING AN INPUT AND AN OUTPUT FILE 

CALL EXCMS (’FILEDEF 9 DISK RESULT DATA Al') 

CALL EXCMS ('FILEDEF 15 DISK FRESULT DATA Al') 



* 



10 



READING ESTIMATED RELIABILITIES 
DO 10 1=1,5000 
READ(9,*) A(I) 

CONTINUE 



SORTING OPERATION 
N=I-1 
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DO 30 I=N,2,-1 
FLAG='OFF' 

DO 20 J=1,I-1 

IF (A(J).GT.A(J+D) THEN 
TEMP=A(J) 

A(J)=A(J+1) 

A(J+1)=TEMP 

FLAG='0N' 

END IF 

20 CONTINUE 

IF ( FLAG. EQ. 'OFF') THEN 
GO TO 40 
END IF 

30 CONTINUE 
40 CONTINUE 

#V >** Vv #\i #V iV A A ^ ^ iV A «'» iV iir* v'% A A ?V v*% iV /'* A vV /» V* V» V* j-V «V iW* A A A A v*\ A ^ J% <V A 

WRITING THE RESULTS IN ASCENDING ORDER TO OUTPUT FILE 
DO 50 1=1,5000 

WRITE (15,*) A(I) 

50 CONTINUE 

WRITE (15,1) A(250) 

1 FORMAT (///,15X,'95 LOWER CONFIDENCE BOUND IS ' , 1X,F12 . 7) 

STOP 

END 
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APPENDIX D. PROGRAM INITIAL 



PROGRAM INITIAL 



THIS PROGRAM, CALCULATES INITIAL GUESSES FOR PARAMETERS IN 
LOGLINEAR MODEL BY MEANS OF PROGRAM PARAM WHICH IS IN APPENDIX E. 
PROGRAM SUPPLIES PARTIAL SUMS OF EXPECTATION, TO SOLVE EQUATIONS 
IN PROGRAM PARAM. IT READS INTERACTIVELY NUMBER OF FAILURES IN 
TESTS. PROGRAM WRITES RESULTS TO AN OUTPUT FILE CALLED EXPECT. 



TL*J _f f f f f * f f f f_ < f - ^ - - - -t - - 



VARIABLES 



FOM 

FOTH 

FOV 

FOA 

PI 

P2 

P3 

P4 

Q1 

Q2 

Q3 

Q4 

XI 

X2 

X3 

X4 



NUMBER OF FAILURES IN MANUFACTURER TEST 

NUMBER OF FAILURES IN TEMPERATURE AND HUMIDITY TEST 

NUMBER OF FAILURES IN VIBRATION TEST 

NUMBER OF FAILURES IN ALTITUDE TEST 

SUCCESS RATIO FOR MANUFACTURER TEST 

SUCCESS RATIO FOR TEMPERATURE AND HUMIDITY TEST 

SUCCESS RATIO FOR VIBRATION TEST 

SUCCESS RATIO FOR ALTITUDE TEST 

FAILURE RATIO FOR MANUFACTURER TEST 

FAILURE RATIO FOR TEMPERATURE AND HUMIDITY TEST 

FAILURE RATIO FOR VIBRATION TEST 

FAILURE RATIO FOR ALTITUDE TEST 

PARTIAL SUM OF EXPECTED NUMBERS IN CELLS FOR EQUATION 1 

PARTIAL SUM OF EXPECTED NUMBERS IN CELLS FOR EQUATION 2 

PARTIAL SUM OF EXPECTED NUMBERS IN CELLS FOR EQUATION 3 

PARTIAL SUM OF EXPECTED NUMBERS IN CELLS FOR EQUATION 4 



TYPE DECLARATION 
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REAL F0M,F0TH,F0V,F0A,P1,P2,P3,P4,Q1,Q2,Q3,Q4,PR0D,X1,X2,X3,X4 









OPENING AN OUTPUT FILE 

CALL EXCMS (’FILEDEF 13 DISK EXPECT DATA Al') 



* 

* 



* READING NUMBER OF FAILURES IN EACH TEST INTERACTIVELY 
WRITE(*,*) 'PLEASE WRITE NUMBER OF FAILURES IN MANUFACTURER TEST' 
READ(*,*) FOM 

WRITE(*,*) 'PLEASE WRITE NUMBER OF FAILURES IN TEMP. AND HUM. TEST' 
READ(*,*) FOTH 

WRITE(*,*) 'PLEASE WRITE NUMBER OF FAILURES IN VIBRATION TEST' 
READ(*,*) FOV 

WRITE(*,*) 'PLEASE WRITE NUMBER OF FAILURES IN ALTITUDE TEST' 
READ(*,*) FOA 

* CALACULATION OF SUCCESS RATIOS 
P1=(20.0-FOM)/8.00 
P2=(20.0-FOTH)/8.00 
P3=(32.0-FOV)/8.00 
P4=(20.0-f6a)/8.00 

CALCULATION OF FAILURE RATIOS 
Ql=F0M/8.00 
Q2=FOTH/8.00 
Q3=FOV/8 . 00 
Q4=FOA/8 . 00 

* CALCULATION PARTIAL SUMS OF EXPECTATION FOR PROGRAM PARAM 
PR0D=2 . 0*P3+2 . 0*Q3+2 . 0*P4+2 . 0*Q4 

Xl=4 . 0*Pl+4 . 0*P2+PROD 
X2=4 . 0*Pl+4 . 0*Q2+PR0D 
X3=4 . 0*Ql+4 . 0*P2+PR0D 
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X4=4 . 0*Q 1+4 . 0*Q2+PROD 



WRITING THE RESULTS TO AN OUTPUT FILE 

WRITE (13,*) X1,X2,X3,X4 

STOP 

END 
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